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ABSTRACT 


An extensive examination of NACA Report No. 496 (NACA 496), "General Theory of Aerodynamic 
Instability and the Mechanism of Flutter," by Theodore Theodorsen, is described. The examination 
included checking equations and solution methods and re-computing interim quantities and all 
numerical examples in NACA 496. The checks revealed that NACA 496 contains computational shortcuts 
(time- and effort-saving devices for engineers of the time) and clever artifices (employed in its solution 
methods), but, unfortunately, also contains numerous tripping points (aspects of NACA 496 that have 
the potential to cause confusion) and some errors. The re-computations were performed employing the 
methods and procedures described in NACA 496, but using modern computational tools. With some 
exceptions, the magnitudes and trends of the original results were in fair-to-very-good agreement with 
the re-computed results. The exceptions included what are speculated to be computational errors in 
the original in some instances and transcription errors in the original in others. Independent flutter 
calculations were performed and, in all cases, including those where the original and re-computed 
results differed significantly, were in excellent agreement with the re-computed results. Appendix A 
contains NACA 496; Appendix B contains a Matlab ® program that performs the re-computation of 
results; Appendix C presents three alternate solution methods, with examples, for the two-degree-of- 
freedom solution method of NACA 496; Appendix D contains the three-degree-of-freedom solution 
method (outlined in NACA 496 but never implemented), with examples. 


I. INTRODUCTION 

In a year 2000 engineering note, Zeiler [1] pointed out that several of the foundational papers and texts 
in aeroelastic flutter [2 through 6] contain numerical errors in some of their example problems. It is not 
surprising that such errors exist because - especially in the cases of references 2, 3, and 4, written in the 
1930's and 40's - all calculations were performed "by hand" with pencil and paper, slide rules, and, after 
they were invented, large table-top four-function electro-mechanical calculators. 

Because these foundational papers and texts are often used in graduate courses on aeroelasticity, Zeiler 
recommended that an effort be undertaken to employ the computing resources available today to re- 
calculate the example problems and to publish the results so as to provide a complete and error-free set 
of example problems. Following Zeiler's recommendation, the purpose of this paper is to provide re- 
computed results for the example problems presented in reference 2: NACA Report No. 496, "General 
Theory of Aerodynamic Instability and the Mechanism of Flutter," by Theodore Theodorsen (referred to 
hereinafter simply as "NACA 496"). Future papers will provide re-computed results for references 3 and 
4. 

In order to perform the re-computation of its numerical examples, an extensive study of NACA 496 was 
undertaken by the present author. The study entailed many parts, including reading and re-reading the 
paper, checking and re-checking its equations, and checking and re-checking its solution method. In the 
process, the present author discovered in NACA 496 computational shortcuts (time- and effort-saving 
devices for engineers of the time) and clever artifices employed in its solution methods. Unfortunately, 
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the present author also discovered, via tripping, several tripping points (aspects of NACA 496 that have 
the potential to cause confusion) and errors. These items are discussed in the present paper. 

In reference 1, Zeiler stated "One does not set about lightly to correct the masters ..." Embracing this 
notion, the present results have been carefully checked and re-checked, and, as a further check, for 
some of the numerical examples, independent calculations were performed by a colleague. In all cases, 
including those where the original and present results differed significantly, the independent 
calculations were in excellent agreement with the present results, providing additional confidence in the 
present results. 

The arrangement of the remainder of the present paper is as follows: 

Section II contains preliminaries, which are intended to aid in the understanding of the 
remainder of this paper. 

Section III contains an overview of NACA 496, including its theoretical development and solution 
method. 

Section IV describes the checks performed by the present author on the equations and solution 
method of NACA 496. 

Section V presents the re-computation of quantities contained in Tables l-IV of NACA 496. 

Section VI presents the re-computation of the example problems contained in Appendix II of 
NACA 496. 

Section VII contains concluding remarks. 

Appendix A contains a reproduction of NACA 496. (Throughout the present paper the reader 
will be referred to Appendix A frequently.) 

Appendix B contains a Matlab® program that implements the two-degree-of-freedom flutter 
equations and the two-degree-of-freedom solution method found in Appendix I of NACA 496. 

Appendix C contains descriptions of, and example problems for, three alternate solution 
methods for the two-degree-of-freedom flutter equations presented in NACA 496. 

Appendix D contains a description of, and example problems for, the three-degree-of-freedom 
flutter equations presented in NACA 496. 
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II. PRELIMINARIES 


As stated in the Introduction of the present paper, the information contained in this section is intended 
to aid in the understanding of the remainder of this paper. 


Terminology 

The present paper refers to sections and appendices within itself and sections and appendices within 
NACA 496. The present paper also refers to the authors of both papers. To avoid confusion in regard to 
which paper or author is being referred to, the following (unfortunately but necessarily) awkward- 
sounding convention is adopted: for the present paper, the terms "present paper," "this paper," "main 
body," "Appendix A, B, C, or D of the present paper," and "present author" are used; for NACA 496, the 
terms "NACA 496," "main text," "Appendix I, II, or III of NACA 496," and "author of NACA 496," are used. 


Symbology 

The present paper retains the symbology of NACA 496. Appendix A of the present paper (the 
reproduction of NACA 496) contains a symbol list (pp. 9 and 10 of NACA 496). Symbols not contained in 
the list are defined as they are used in NACA 496 and as they are used in the present paper. 


Tripping Points and Errors 

As used in the present paper, the term "tripping point" is defined to be an aspect of NACA 496 that has 
the potential to cause confusion for a reader of NACA 496. Tripping points that the present author 
encountered while reading NACA 496 are identified throughout the present paper and are collected in 
Table I of the present paper. 

Errors in some of the equations, tables, and figures of NACA 496 were found by the present author. 
These errors are identified throughout the present paper and are collected in Table II of the present 
paper. 


Three versions of NACA 496 

The author of the present paper had access to three different versions of NACA 496. At the conclusion 
of their respective main texts, all three versions contain the following lines of text 

"Langley Memorial Aeronautical Laboratory, 

National Advisory Committee for Aeronautics, 

Langley Field, Va., May 2, 1934." 
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indicating that all three are "the same" report. This common date is important, as will be seen below. 

The first version of NACA 496 was in hard-copy form and was contained in the bound volume 
"Twentieth Annual Report of the National Advisory Committee for Aeronautics," which contains the 
significant findings of NACA for fiscal year 1934 (July 1933 through June 1934). The page numbering 
(pages 413 to 433) in this version reflects the fact that NACA 496 was one of many papers in the bound 
volume. Errors in some of the equations, tables, and figures were found (by the present author) to exist 
in this version. 

The second version was in electronic form and was acquired from the NASA Technical Report Server 
website (http://ntrs.nasa.gov). This version was a reprint of the first version and was issued as a stand- 
alone report. The date on the cover of this version is 1949. The page numbering in this version (pages 3 
to 23) reflects the fact that it was a stand-alone report. Some of the errors found in the first version are 
corrected in the second version. 

The third version was also a hard-copy version, was also a stand-alone report, and was acquired from 
the NASA-Langley technical library. Pages 3 to 23 in this version are identical to those in the second 
version, but the date on the cover of the third version is 1940. 

The first tripping point for the present author was the fact that there are three versions of the paper, 
each with a different publication date. A second tripping point was that some of the errors in the first 
version were corrected in the second and third versions, but no known errata was issued by NACA, as 
confirmed by the staffs of the NASA-Langley and NASA-Ames technical libraries. (Errors from the first 
version that were not corrected in the second and third versions, and therefore remain in the second 
and third versions, will be discussed throughout the present paper.) 

The third version, a corrected version, is contained in Appendix A of the present paper. 


Referring to Equations 

NACA 496 contains many equations, but only a minority is identified by letters or numbers. In the main 
text, capital Roman letters and Roman and Arabic numerals, all within parentheses, are used as 
identifiers. In Appendix I of NACA 496, Arabic numerals within parentheses are used as identifiers; in 
addition, some terms are identified by Arabic numerals, but without parentheses. In what would be 
very unconventional by today's report-writing standards, within NACA 496, there are multiple instances 
of different equations being identified using the same equation number. 

In the present paper, references to specific equations in NACA 496 will cite the NACA 496 page and 
equation numbers. For example, equation (A) on page 10, equation (XXI) on page 12, and equation (9) 
on page 14 would be identified herein as "equation (10/A)," "equation (12/XXI)," and "equation (14/9)," 
respectively. 

Equations in the present paper are identified by Arabic numerals within parentheses. 
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Computational Shortcuts 


Given the very rudimentary computational aids of the early 1930's (paper and pencil, slide rules, four- 
function electro-mechanical calculators), it is natural to deduce that in the construction of their solution 
methods engineers of the time would have sought ways to minimize the overall number and/or to 
simplify the complexity of calculations so as to minimize the human time and effort required to obtain a 
solution. It should not be surprising, then, that such time- and effort-saving techniques are evident in 
the NACA 496 solution method. These techniques will be given the general term "computational 
shortcuts." They will be identified throughout the present paper and are collected in Table III of the 
present paper. 


III. OVERVIEW OF NACA REPORT No. 496 

NACA 496 lays out the theoretical development of aeroelastic flutter for a typical section with degrees 
of freedom in torsion (a), aileron deflection (jB), and vertical deflection (at times referred to in NACA 496 
as flexure) (/?). NACA 496 has three appendices: Appendix I presents a detailed description of the steps 
of the implementation of the solution procedure, including tables of numerical values of quantities 
required in the solutions; Appendix II presents a number of numerical calculations, including a few 
comparisons with experimental data; and Appendix III presents integral evaluations for some of the 
velocity potentials. 


Theoretical Development 

In NACA 496, the theoretical development begins with four simplifying assumptions: (1) the flow is 
potential and non-stationary; (2) the "wing" is actually a two-dimensional typical section with no 
thickness and therefore with no airfoil shape; (3) the "wing" motions are sinusoidal and infinitesimal; 
and (4) the "wing" has no internal or solid friction, resulting in no internal damping forces. 

Next, the circulatory and non-circulatory velocity potentials are developed. The magnitude of the 
circulation, in the form of Theodorsen's circulation function, is developed. After this, the aerodynamic 
forces and moments are determined via the chordwise integration of the velocity potentials. Then, the 
aerodynamic forces and moments are combined with the inertia and restraining forces and moments to 
produce three second-order differential equations in the three unknowns a, /?, and h, equations (10/A), 
(10/B), and (10/C), reproduced here: 
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where C(k), Theodorsen's circulation function, is a complex function of reduced frequency, k. Equation 
(10/A) defines the sum of the moments about the elastic axis; equation (10/B), the sum of the moments 
about the aileron hinge; and equation (10/C), the sum of the forces on the entire "wing" in the vertical 
direction. Upon examination of equations (10/A) through (10/C) one sees that k appears only implicitly 
gin C(k) while velocity, v, appears multiple times. 

Figure 1 in the present paper is a re-drawn figure 2 from NACA 496. It illustrates the definitions and 
positive senses of many of the important physical parameters appearing in the equations of motion. 


Assumed forms of the unknowns, a, B, and h . - The assumed forms of the unknowns in equations 
(10/A) through (10/C) are introduced in NACA 496 as sine functions of the distance, s, traveled by the 
wing from the first vortex element. Employing the distance formula, s = vt, the sine functions in 
complex form may be expressed as: 


a 

II 

a 

o 

sr 
cy| ^ 

Ct 

(la) 

p = 

p 0 e lik ^ +<Pl) 

(lb) 

h = 

/i 0 e i(fe F t+< P 2 ) 

(lc) 


where a 0 , f3 0t and h 0 are the (infinitesimal) amplitudes of a, /?, and h, and cp ± and cp 2 are phase angles of 
/?and h with respect to a. The first and second time derivatives of a, /?, and h are: 


a = ik-a and a = 

b 


/? - ik-/3 and /? = 


h = ik-h and h = 

b 


~( k i ) “ 

-(* if » 
-i k if h 


Substitution of assumed forms into, and normalization of, equations . - Making the substitutions of 
equations (la) through (lc) and their time derivatives into equations (10/A) through (10/C) transforms 
the latter equations from three simultaneous differential equations into three simultaneous algebraic 
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equations. The algebraic equations are then normalized by the quantity k \j~kj and by their 
respective exponentials and amplitudes, resulting in the equations taking the form 


(Aaa + £l a X) a + AapP + A ah h = 0 

(2a) 

Aba a + (A b p + Q. b X){S + A bh h = 0 

(2b) 

A ca a + A c p/i + ( A ch + £l h X)h = 0 

(2c) 


The physical constants of the problem ( k, q, b, c, Xa, xp f rp ; coa, cop, and coh) reside in the coefficients, 
that is in the A's and the fix terms. 

The quantities A aa , etc. are complex functions of reduced frequency, k, with real parts R oa , etc. and 
imaginary parts l aa , etc. The real parts are defined in equations (12/1) through (12/9) and the imaginary 
parts in equations (12/11) through (12/19). (There is no equation (12/10). This equation numbering 
convention was chosen in NACA 496 so that for any given imaginary part, the second digit in its two-digit 
equation number would be the same as the single-digit equation number of its corresponding real part.) 

The QX terms are discussed next. 


Products f2aX, f2nX, and fihX . - In equations (2a), (2b), and (2c) the products fl a X, flpX, and I2hX are 
real. They are derived from quantities (from eqn. (10/A)), (from eqn. (10/B)), and ^(from 
eqn. (10/C)), respectively. Via the substitution and rearrangement of terms and the use of cancelling 
expressions in the numerator and denominator, from page 12 of NACA 496 these products are 


r 

_ 

/"aV 

\ 2 h 

(br r co r \ 2 

(3a) 

k 2 Mv 2 K 

\ (o r r r , 

) K 1 

\ vk / 

Cp _ 

fuprp' 

f-l 
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(3b) 

k 2 Mv 2 K 

\ a > r r r , 

1 K 

\ vk ) 

C h b 2 

( "ft > 

| 2 -( 

'br r (o r \ 2 

(3c) 

k 2 Mv 2 K 

\(o T r T ) 

1 k' 

\ vk ) 


where, to the right of the second equal sign in each equation, X comprises the two right-most terms 


'■Kt) 

and the respective Q’s, the remaining terms 



(4) 


(5a) 
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(5b) 



The quantities co r and r r are a reference frequency and reference length, which may be conveniently 
chosen. 


Referring back to the normalization that produced equations (2a) through (2c), a critically important 
result of the normalization is that the quantity X has been conveniently isolated from the other terms in 
these equations and X is the only quantity in these equations that contains the velocity. The NACA 496 
solution methods take advantage of this result in an ingenious way. As will be seen immediately below, 
X and the Q ' s are clever artifices created by the author of NACA 496 to enable the determination of the 
flutter velocity, v/, and flutter reduced frequency, kf. At times X and the /2's are treated not as the 
known quantities expressed in equations (4) and (5), but rather as parameters. At other times X and the 
/2's are treated as the known quantities. 


Solution Methods 


The solution of equations (2a), (2b), and (2c) is obtained when their determinant is zero 


A aa + ft a X A ap 

A bcc A b(l + 

A ca 


A ah 

A bh 

A ch + 


= 0 


( 6 ) 


The form of equation (6), det(A + XQ) = 0, bears some resemblance to the characteristic polynomial of 
the eigenvalue problem, det(A - XI) = 0. The quantity X in eqn. (6) would be analogous to X in the 
eigenvalue problem and the quantity Q, a 3x3 diagonal matrix with diagonal elements Q a , £2p, and f2 h , 
would be analogous to the identity matrix. However, in terms of the solution of equation (6), this is 
where its resemblance to the eigenvalue problem ends. 

Expanding the determinant in equation (6) yields equation (12/XXI), re-written here 

"1" "f" "f" ^h^aAb 

+(n a M aa + + n h M ch )X + D = 0 

which is third-order in X with complex coefficients. In equation (7), D is the determinant in equation (6) 
without the QX terms and the M's are the minors of the determinant without the fXX terms. Except for 
Q a , Qp, and Qh, all quantities in equation (7) are functions of reduced frequency, k. 
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The object of the solution method is to find the values of v and k that cause equation (7) to be satisfied. 
Such values are then identified as v/ and kf, the flutter velocity and the flutter reduced frequency, 
respectively. 


Three-degree-of-freedom solution method. - NACA 496 outlines the method by which the three- 
degree-of freedom flutter problem (that it names the "general case") may be solved. 

NACA 496 does not recommend solving directly the third-order equation with complex coefficients (eqn. 
(7)). Instead, it suggests the computational shortcut of creating two equations, each with real 
coefficients, by separating the real and imaginary parts of the original. The first equation is cubic in X 
and is obtained from the real parts of the A's, M's, and D; the second is quadratic in X and is obtained 
from the imaginary parts. The right hand side of both equations is zero. This separation of a complex 
equation into its real and imaginary parts eliminates complex arithmetic from the solution method. 


At this point in the solution method, the artifice of treating X as a parameter, rather than as a known 
quantity, is employed. The cubic and quadratic equations are each solved for their roots, identified 
herein as X Ri ,X R2 ,X R3 ,X Iif and Xj 2 for a large number of reduced frequencies. The X R 's and X's are 
then plotted on the same set of axes as functions of the inverse of reduced frequency. Intersections of 
any of the X R 's with any of the X! s identify the values of X and the values of reduced frequency, which 
simultaneously satisfy the cubic and quadratic equations, and therefore also solves the original third- 
order equation with complex coefficients. These values are the flutter values -Xf and kf. 


With Xf and the /(/known, the artifice of treating X as a parameter is now abandoned and equation (4) is 
employed, which when re-written to solve for flutter velocity, becomes 


1 1 bco r r r 

Vf = 


( 8 ) 


where /cand b are known constants of the problem and co r and r r are chosen reference values. 


Although many illustrative examples are presented in Appendix II of NACA 496 none are given for the 
general case, almost certainly due to the difficulty in the early 1930's of solving "by hand" a cubic 
equation. 


Two-degree-of-freedom solution method. - There are three two-degree-of-freedom subcases available 
from the three-degree-of-freedom equations: torsion - aileron; aileron - flexure; and flexure - torsion. 
NACA 496 lists the subcases in the order just given, but then identifies the first one listed as subcase 3, 
the second one listed as subcase 2, and the third one listed as subcase 1. No reason is offered in NACA 
496 for this curious pairing of listing order and subcase number. 

NACA 496 presents an ingenious solution method for obtaining the flutter velocity and flutter reduced 
frequency for the three two-degree-of-freedom subcases. This solution method again involves 
computational shortcuts and the artifices, but, in yet another curiosity, this solution method differs 
fundamentally from the recommended solution method for three degrees of freedom. 
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Whereas equation (12/XXI) is the characteristic equation obtained from the expansion of the 3x3 
characteristic determinant (eqn. (7) in the present paper), equations (12/XXI I), (12/XXIII), and (12/XXIV) 
are the characteristic equations obtained from the expansion of the three 2x2 minors of the 
characteristic determinant. The reader is referred to Appendix A of the present paper. 

Using the torsion - aileron subcase (subcase 3) as the example, NACA 496 goes through the algebraic 
steps of this solution method. The algebraic steps for subcases 2 and 1 mirror those for subcase 3. 
Equation (12/XXII), a second-order equation in X with complex coefficients, is the starting point. NACA 
496 again employs the computational shortcut of creating two equations, each with real coefficients, by 
separating the real and imaginary parts of the original. The first equation is quadratic in X and is 
obtained from the real parts of A aa , Abp, and /VU; the second is linear in X and is obtained from the 
imaginary parts. The right hand side of both equations is zero. 

At this point the two-degree-of-freedom solution method departs from the three-degree-of-freedom 
solution method, but like the latter solution method it also employs the artifices: it treats both X and Q a 
as parameters, rather than as the known quantities in equations (4) and (5a). The quadratic and linear 
equations and parameters X and Q a are treated as two equations in two unknowns and are solved 
conventionally by substitution. 

From page 13 in NACA 496, the first step is to solve the linear equation forX, yielding 


+ ttplaa 

This expression is then substituted into the quadratic equation in X, thereby eliminating X from that 
equation and producing, after several steps, equation (13/XXV), which is quadratic in i2«and re-written 
here 


^ai^^ch^bp M C hRbpIbp > ) + n a [ M C b(R a pIba + lafiRba) + 2M c ^/ aa /^] 

“h ^ch^aa ~ ^ch^aa^aa ~ 0 

Equation (10) is then solved for its two roots, H ai and ft a2 , for a large number of reduced frequencies. 

With the artifice still in place, but recognizing from equation (5a) that Q a is defined as the square of a 
quantity, only the real positive values of H ai and fl a2 and their corresponding reduced frequencies are 
retained. A plot of Q a as a function of 1/k, which contains both fl^and Cl a2 , is produced. 

Next, values of X are found by substituting the values of Q a , just obtained, their corresponding values of 
reduced frequency, and the values of M' C h, laa, and I bp, at corresponding reduced frequencies, into 
equation (9). NACA 496 chooses to present the quantity X as a function of Q a rather than as a function 
of 1/k. 

At this point a new quantity, F, the nondimensional flutter factor, is introduced 
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1 

X 


( 11 ) 


F 


1 

k 


F is also treated as a parameter. Substituting into equation (11) the X's and their corresponding reduced 
frequencies produces a plot of Fas a function of Q a , again containing both fl^and O a2 . 


The two curves, C2 a vs. 1/k, and Fvs. Q a , represent a family of flutter solutions for a range of C2 a 's. Each 
point on the first curve has a corresponding point on the second and each pair of corresponding points 
represents a unique flutter solution. 


To obtain a specific flutter solution from the family of solutions, the artifice is abandoned: quantities 
Q a , X, and F are no longer interpreted as parameters and equations (5a) and (4) are employed. In 
addition, when the expression in equation (4) is substituted for X in equation (11), the flutter factor 
becomes 


\[kv 

ba) r r r 

When solved for flutter velocity, equation (12) becomes 


v f 


= F 


ba) r r r 


( 12 ) 


(13) 


Thus, once a family of flutter solutions is obtained, the following simple two-step process for obtaining a 
specific flutter solution from the family of solutions is performed: 


First, problem-specific values of co a , r a , a> r , and r r are substituted into equation (5a), yielding 
a problem-specific value of Q a - This value of Q a is located in the plot of Q a vs. 1/k. The 
value of 1/k corresponding to this value of Q a then yields, via its inverse, the flutter reduced 
frequency, kf. 


Second, the problem-specific value of Q a is located in the plot of Fvs. Q a - The value of F 
corresponding to this value of Q a and problem-specific values of b, co r , r r , and /care 
substituted into equation (13), yielding the flutter velocity, v/. 


IV. CHECKS OF NACA 496 EQUATIONS AND SOLUTION METHOD 

Beginning with equations (10/A), (10/B), and (10/C) in the "DIFFERENTIAL EQUATIONS OF MOTION" 
section of NACA 496, the present author checked all of the equations in the main text and Appendix I. 
With a few exceptions, discussed below, all equations checked. 
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Missing 1/b Factors in Expressions on Page 12 

In checking the equations in NACA 496, the present author duplicated 12 of the 18 expressions in the 
left column on page 12. However, as they appear in NACA 496, the real expressions in equations (12/3), 
(12/6) and (12/9) and the imaginary expressions in equations (12/13), (12/16) and (12/19) each lack a 
factor of 1/b compared to what the present author found. These six expressions make up the entirety 
of equation (10/C) and none of the six appears in either equation (10/A) or equation (10/B). 

As it turns out, the missing 1/b factors are of no consequence because the left-hand side of an equation 
whose right-hand side is zero, such as equation (10/C), may be multiplied by any factor without 
changing the equation. The present author believes that, as a computational shortcut, the author of 
NACA 496 multiplied both sides of equation (10/C) by the semi-chord, b, thus eliminating the 1/b factor 
from equations (12/3), (12/6), (12/9), (12/13), (12/16), and (12/19) and thereby relieving engineers of 
the time of extra, but unnecessary, multiplications in the hand calculations of these expressions. 


Comparing Real Expressions in Main Text and Appendix I 

The real expressions in the left column on page 12 (main text), the real expressions in the right column 
on page 14 (Appendix I), and the real expressions in the left column on page 15 (Appendix I), term-by- 
term, share the same equation numbers. This assignment of equation numbers is a potential tripping 
point because, term-by-term, the expressions obviously differ from each other. As it turns out, the 
expressions on page 12 are, term-by-term, equal to the sums of the respective expressions on pages 14 
and 15. In checking the equations the present author reproduced all of these expressions. 

Table III on page 17 in Appendix I contains a range of values for the expressions in equations (14/1) 
through (14/9). 


Comparing Imaginary Expressions in Main Text and Appendix I 

The imaginary expressions in the left column on page 12 (main text) correspond to the imaginary 
expressions in the right column on page 14 (Appendix I). A cursory examination of these equations 
suggests the identity of corresponding expressions because: (1) the same symbols for the I's (l aa , kp, Ich, 
etc.) are used in both places; and (2) the same equation numbers are used in both places. 

In checking the equations, the present author reproduced the expressions on page 12. However, in 
comparing the expressions on page 12 with those on page 14, it is seen that each expression on page 12 
contains a multiplying factor of 1/k while those on page 14 do not. As with the omission of the factor of 
1/b previously discussed, it turns out that the omission of 1/k is also ultimately of no consequence, and 
for the same reason. As with the 1/b omission, the present author believes that the 1/k omission was 
also done deliberately by the author of NACA 496 as another computational shortcut, again relieving 
engineers of an extra, but unnecessary, multiplication in the hand calculation of these expressions. 
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The author of NACA 496 does alert the reader to these differences in the I's, but via a whisper, not a 
shout: a footnote on the bottom of page 14. No explanation is given as to why 1/k has been eliminated, 
yet the elimination has important consequences in the equations employed in the solution method. The 
differences represented a significant tripping point for the present author, who initially missed the 
footnote. 

Table IV on page 17 in Appendix I contains a range of values for the expressions in equations (14/11) 
through (14/19) (that is, the expressions without the 1/k factor). 


Comparing Equations for Solution Method in Main Text and Appendix I 

The NACA 496 solution method is derived beginning at the bottom of the left column on page 12 and 
continuing through the end of the main text on page 13. The procedure for implementing the solution 
method is contained in Appendix I on pages 14 through 16. These parts of NACA 496 will be referred to 
as the derivation section and the implementation section, respectively. 

The solution method involves solving quadratic equations whose coefficients (a, b, and c) contain 
numerous sums and differences of products of the R's and I's and other terms. The other terms are the 
real and imaginary parts of the minors of the characteristic determinant (M R C h, M R aa , M R bp, and M' C h, 

M'aa, M'bp, and referred to hereinafter as "real minors" and "imaginary minors"), which also contain 
sums and differences of products of the R's and I's. 

As mentioned above, to take advantage of a computational shortcut, all the I's in the implementation 
section had their common 1/k factor eliminated. For clarity these I's without the 1/k factor will be called 
"new I's." The new I's and any quantities involving the new I's (the real and imaginary minors and the 
quadratic coefficients a, b, and c) are all affected by the elimination of the 1/k factor. 

This situation would have caused no confusion had new symbols (perhaps original symbols with tildes) 
been introduced in NACA 496 to represent the new I's. The original symbols could have been used in 
the derivation section; the new symbols in the implementation section. However, the author of NACA 
496 chose not to do this; he chose to retain the original symbols in the implementation section and in so 
doing assigned each symbol two meanings, with each meaning specific to a given section of the paper. 

However, with careful attention to section-specific definitions and by doing a little algebra, it is easy to 
prove that the quadratic coefficients in the implementation section are proportional, term-by-term, to 
the quadratic coefficients in the derivation section, with proportionality constant 1/k 2 . Thus, because 
multiplying either equation (both of whose right-hand sides are zero) by the proper quantity (either k 2 
or 1/k 2 ) will produce the other, both equations will yield the same roots, and therefore the same values 
of v/ and kf. 
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V. RE-COMPUTATION OF QUANTITIES EMPLOYED IN NACA 496 SOLUTION METHOD 


NACA-496 contains four tables of quantities needed in the implementation of the NACA-496 solution 
method. For the engineer of the time, these tables, especially the latter two, would have saved valuable 
human computation time. 


Constants Resulting from the Integration of Velocity Potentials 

Table I in NACA 496 contains constants, 7„ associated with integrations of the velocity potentials. These 
constants are functions of the nondimensional chordwise location of the aileron hinge, c. At the top of 
the right column on page 5 of NACA 496 are the integrals in question and it can be seen that there are 
two types of integrals that produce the 7/ s: those whose limits of integration are c to 1 (integrations 
over the aileron only) and those whose limits are -1 to 1 (integrations over the entire "wing"). 

At the bottom of the right column on the same page are the equations for Ti through Ti 4 . However 
Table I lists only 11 of the 14 7/s. The omitted 7/s (Tg, T 13 , and T 14 ) are also functions of a, the elastic axis 
location. 

Figure 2 in the present paper contains comparison plots of the 7/s listed in Table I and present 
calculations of the 7/s. The open circles correspond to the tabular values, which are shown in NACA 496 
for five values of c (-1, - Yi, 0, 1)- The solid lines correspond to the present calculations, performed in 

Matlab® and computed from c = -1 to c = 1 in increments of 0.01. The plots show there is excellent 
agreement between the original results and the present results. 


Theodorsen Circulation Function 

Table II in NACA 496 contains the real part, F, and the negative of the imaginary part, -G, of the 
Theodorsen circulation function and their component Bessel functions. All quantities are functions of 
reduced frequency, k, and are shown for 16 values of k. 

Figure 3 in the present paper contains comparison plots of F and -G listed in Table II and present 
calculations of 7 and -G. The horizontal axis is logarithmic to more clearly show (compared to what a 
linear axis would show) the comparisons at low values of 1/k. A logarithmic axis "spreads out" the 
values of F and -G at the low values of 1/k and compresses them at the high values of 1/k in such a way 
that the distance on the page between 1/k of 0.1 and 0.2 is the same as the distance on the page 
between 1/k of 10 and 20. The open circles correspond to the tabular values; the solid lines correspond 
to the present results, computed in Matlab® from 1/k = 0.1 to 1/k = 40 in increments of 0.004. 

For values of 1/k of 5 and below there is excellent agreement between the original results and the 
present results. For all values of 1/k except 0.1, the original results and the present results differ by less 
than one percent. For values of 1/k of 10 and above the differences vary from 0.2 to 5.6 percent, with 
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the largest differences occurring at 10 and 40. These largest differences will produce corresponding 
differences seen in the illustrative examples. 


Reduced-Frequency-Dependent Portions of the Elements of the Equations of Motion 

Table III in NACA 496 contains the frequency-dependent portions of the real parts of the coefficients of 
the equations of motion; Table IV contains the imaginary parts. All quantities in Tables III and IV are 
functions of reduced frequency, k. Some quantities in Tables III and IV are functions of the 
nondimensional chordwise location of the elastic axis, a; other quantities are functions of the 
nondimensional chordwise location of the aileron hinge, c; some quantities are functions of both; and 
two quantities are functions of neither. The quantities in Tables III and IV will be discussed together. 

Figure 4 in the present paper contains comparison plots of the real and imaginary parts from Tables III 
and IV and from present calculations, again performed in Matlab®. As was the case for figure 3, in this 
figure the horizontal axis is logarithmic to more clearly show (compared to what a linear axis would 
show) the comparisons at low values of 1/k. Figure 4(a) contains elements from the pitching moment 
equation; figure 4(b), from the hinge moment equation; and figure 4(c), from the vertical force 
equation. For the present calculations the inverse of reduced frequency ranged from 1/k - 0.1 to 1/k - 
40, in increments of 0.0004. 

The plots are for values of a and c of -0.4 and 0.5, respectively. Except for the comparison plot for R'^ h , 
which will be discussed in the next paragraphs, all comparisons show reasonably good agreement. The 
comparisons at values of 1/k of 10, 20, and 40 consistently show differences and are due to the 
differences in Fand-G, discussed above. 

The comparison plot for R'^ h in figure 4(a) shows significant differences between the original and 
present results. Close examination of the values in Table III reveals that, for a = -0.4, R'^ h is very nearly 
directly proportional to 1/k, that is, is very nearly linear with 1/k. To more easily see this linearity, 

R'a h as a function of 1/k has been redrawn in figure 4(d) with a linear horizontal axis. A red dashed line 
has been added to emphasize the linear relationship between R^ h and 1/k. 

If one assumes the values of R'^ h in Table III are correct and then employs equation (14/3) to back 
calculate to find the values of G necessary to produce these values of R'^ h , one finds that the resulting 
G's are independent of 1/k and equal to the value of G at 1/k = 1. (According to Table II in NACA 496 
and figure 3 in the present paper, G varies from about -0.2 to 0 and clearly is a function of 1/k.) The R’^ h 
row of Table III corresponding to a = -0.4 and c = 0.5 is in error. 


VI. RE-COMPUTATION OF ILLUSTRATIVE EXAMPLES IN NACA 496 

This section of the paper presents the re-computation of the illustrative examples contained in Appendix 
II of NACA 496. These examples comprise all three two-degree-of-freedom subcases that may be 
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derived from the three-degree-of-freedom general case. NACA 496 presents no examples for the 
general case. 

A brief orientation to figures 5 through 17 in the present paper: figure-by-figure they contain scanned 
reproductions of the original figures 5 through 17 from NACA 496, but with the addition of the present 
re-computations. In all cases, the original results are in black and the re-computed results are in color. 
Because of "fuzziness" resulting from the scans, the figure legends for these figures have been re-typed. 

Additionally, in figures 5 through 11 and 15 through 17, the colored O and + symbols are results from 
the present implementation, but using the tabular values in Tables I, II, III, and IV of NACA 496, and the 
colored solid and dashed lines are results from the present implementation using present computed 
values throughout. 

As mentioned in Section III of the present paper, the plots of /2vs l//cand Fvs Q are comprised of two 
values of Q, the two quadratic roots obtained from the solution of equation (13/XXV) (as well as two 
values each from the solutions of eqns. (13/XXVI) and (13/XXVII)). With this in mind, the colored solid 
lines and the colored O symbols correspond to the first quadratic root; the colored dashed lines and the 
colored + symbols correspond to the second quadratic root. 


Implementation of NACA 496 Solution Method 

The implementation of the NACA 496 solution method was performed by the present author via the 
writing and execution of m-files in Matlab®. The equations outlined in the main text and Appendix I of 
NACA 496 were solved for 100,001 values of the inverse of reduced frequency, beginning at zero, with 
an increment of 0.0005. 


Standard Case 

Figures 5 through 10 present what NACA 496 calls a "standard" case, represented by the quantities: 

k- 0.1; c = 0.5; a = -0.4; x a = 0.2; r a 2 = 0.25; xp-V^ 0 ; r/ = y l 6 o; and co a , cop, a>h variable. 

Figures 5 and 6 present results for Subcase 3, torsion-aileron; figures 7 and 8 present results for Subcase 
2, aileron-flexure; and figures 9 and 10 present results for Subcase 1, flexure-torsion. 

For these three subcases, the first figure in each pair shows its respective plot of /^versus 1/k; the 
second figure in each pair shows its respective plot of F versus Q. However, the quantity F shown in 
these figures is not the quantity F defined on the bottom of page 15 and the top of page 16 of NACA 496 
(and by eqn. (11) in the present paper). The quantity F shown in figures 6, 8, and 10 is actually the 
square of the defined quantity. This double use of a symbol to represent not only its original definition, 
but also the square of its original definition, represented a significant tripping point for the present 
author. 
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Subcase 3 . - Figures 5 and 6 contain the results for this subcase. While there are some differences 
between the original and present results in these figures, the overall magnitudes and shapes of the 
respective curves are very similar. In figure 5, the agreement between original and present results is 
excellent at values of 1/k greater than 0.5, but progressively less so at lower values of 1/k. In figure 6, 
the agreement is better at the lower values of Q a than at the higher values. In both figures, the present 
results represented by the colored symbols and colored lines exhibit excellent agreement with each 
other, reflecting the excellent agreement between the original tabular quantities and the present 
computed quantities, as illustrated in figure 4. 

Subcase 2 . - Figures 7 and 8 contain the results for this subcase. In each figure, there are three sets of 
results, labeled (a), (b), and (c), representing three different values of xp, the center of gravity of the 
aileron. The original and present curves (a) and (b) in figure 7 agree fairly well in the overall magnitudes 
and shapes of the respective curves. However, the original and present curves (c) do not agree well in 
either magnitude or shape. The original and present curves in figure 8 echo the agreement and 
disagreement noted for the curves in figure 7. In both figures, the present results represented by the 
colored symbols and colored lines exhibit excellent agreement with each other, reflecting again the 
excellent agreement between the original tabular quantities and the present computed quantities, as 
illustrated in figure 4. 

Subcase 1 . - Figures 9 and 10 contain the results for this subcase. Once again, in these figures the 
original and present results are very similar in terms of the overall magnitudes and shapes of the 
respective curves. In figure 9, the original curve does not "touch" the vertical axis, whereas the present 
curve does (as did all original and present curves in figures 5 and 7). This trend is repeated in figure 10, 
but for the horizontal axis. The present results represented by the colored symbols and colored lines do 
not agree as well in figure 9 and 10 as they did in figures 5 and 6 and figures 7 and 8. These differences 
are due in part to the errors in R'^ h in Table III, and illustrated in figure 4, as previously discussed. 


Parameter Variations for Subcase 1 

Figures 11 through 14 contain plots of flutter factor, F, as defined on pages 15 and 16 of NACA 496, as 
functions of various parameters. Figure 11 presents results for variation in frequency ratio, coi/cor, 
figure 12, for variation in the location of the elastic axis, o; figure 13, for variation in the radius of 
gyration, r«; and figure 14, for variation in the location of the center of gravity, x a . All results in this 
subsection of the paper are for Subcase 1, flexure-torsion. 

For figures 12, 13, and 14, as a check on the present results, independent calculations were performed 
by a colleague of the present author who implemented in Matlab® the equations found in reference 5 
for a typical section, including unsteady aerodynamics defined by the Theodorsen circulation function. 
The code loops on velocity. For each velocity, the calculations are run through an iterative loop until the 
frequencies are converged, followed by an evaluation of damping values to determine if the system is 
stable or unstable. 
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To be able to directly compare the independent flutter results with the present flutter results, the 
independent flutter velocities were converted to flutter factors using equation (12). As will be seen 
below, in all cases, the present results and the independent results are in excellent agreement, differing 
in F by no more than 0.5% and in reduced frequency by no more than 3.7%, thereby simultaneously 
confirming and affording confidence in the present results. 

Variation of Frequency Ratio . - The original figure 11 contains theoretical and experimental results, in 
black. The present results are shown in colored symbols and colored curves. The quantities specific to 
this example are: 


k-V^ o 0 ; a = -0.4; x a =0.2; r a - 0.5; coi/co 2 variable. 

This figure is similar to figure 10, but as already stated, the F in this figure is the as-defined F. In 
addition, in this figure, equation 5(c) has been employed to convert the units of the abscissa from Qh to 
( O 1 /CO 2 (that iS, (Dh/CDa)- 

Across the limited range of frequency ratio chosen for this figure, the original theoretical results and 
both sets of present results are in excellent agreement. The colored symbols in the figure correspond to 
values of 1/k of 10, 20, and 40. The re-computed curve in this figure does not "touch" the horizontal 
axis as did the re-computed curve in figure 10 because the much smaller value of mass ratio, k, for 
figure 11 changes the character of the solution. 

Variation of Elastic Axis Location . - In addition to containing the dependency of F on elastic axis 
location, a, figure 12 also contains the dependency of a quantity D on a. According to the sign 
convention established in figure 2 of NACA 496, increasingly negative values of o corresponds to a 
forward movement of the elastic axis. 

The quantities specific to this example are: 

K=y 4 ; coi/ co 2 - %; x a - 0.2; r a - 0.5; a variable. 

In the present effort to re-compute D as a function of a, another tripping point was encountered: D is 
never defined in NACA 496. However, the present author attempted to deduce the definition of D so 
that a comparison could be made between the original and the re-computed curves. 

In the left column on page 16 of NACA 496, there is a short discussion that ends with the equations (un- 
numbered in NACA 496 but numbered herein) 


V D = VrA- (14) 

2 +<X 

v F = v r F (15) 

where a and F are defined as above, Vf is the flutter velocity, vr is a reference velocity defined to be the 
divergence velocity for the specific case of the elastic axis located at the mid-chord (that occurs when a 
= 0 ), and v D is the divergence velocity for other locations of the elastic axis. 
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Velocity v R is a proportionality constant in equations (14) and (15). In equation (15), v R relates the flutter 
velocity, v F , to the flutter factor, F. The present author speculated that in equation (14) v R might relate 
the divergence velocity, v D , to a "divergence factor," D. If this were so, then according to equation (14), 
quantity D would be equal to the inverse of (V 2 + a). As a test of this hypothesis the inverse of (V 2 + a) 
was plotted as a function of a. The hypothesis failed; the curve did not match curve D in figure 12. 
(Equation (14) is clearly incorrect because, according to the definition of v R , when a is zero, v D should be 
equal to v R , but under this condition equation (14) predicts v D to be twice v R .) 


This miss-match led the present author to re-derive the expression relating v D to D, which resulted not in 
equation (14), but in the following expression 


v D = 


VR yfl + 2 a 


(16) 


When the quantity ^== is plotted as a function of a, the colored dashed curve in figure 12 results. 

Except for an area of slight disagreement (less than 5%) between o = -0.2 and a = -0.4, the original and 
present curves agree with each other very well. This overall agreement led the present author to 
conclude that the original curve D in figure 12 represents the expression the author of NACA 496 knew 
to be correct (eqn. (16)) and that the expression on page 16 of NACA 496 (equivalent to eqn. (14)) is a 
misprint. 


The colored circles in figure 12 represent the re-calculation of curve F. Each circle was obtained by 
selecting a value of a and executing the NACA 496 solution method for Subcase 1 in which a curve 
similar to that in figure 10 is used to find the value of F at the following value of Qh - 


n h — 


\ 2 = r o>i \ 2 = fM 2 (W _ m 2 (M 2 = y 

(O a r a J ^u> 2 r a' V r a ) \6 / V0.5/ 


It is obvious that the re-calculation differs markedly from the original. Not only do the numerical values 
significantly differ from each other, but the trends with increasingly negative values of a are opposite 
each other. The three large X's at values of a of -0.2, -0.3, and -0.4 are from the independent calculation 
and are in excellent agreement with the present calculation. The present author can offer no 
explanation for the significant disagreement between the present calculation and the original. 


Variation of Radius of Gyration . - The original figure 13 illustrates the dependency of F on the radius of 
gyration, r a - for two different values of frequency ratio, coi/co 2 . Curve A corresponds to (Oi/co 2 - curve 
B corresponds to coi/co 2 - 1. Other quantities specific to this example are 


K-V^ f a = -0.4; x a = 0.2; /^variable 


For ease of presentation and in an attempt to preserve the figure numbering convention described at 
the beginning of this section, this figure is divided into "offspring" figures: 13A and 13B. Further, figure 
13B has parts (a), (b), and (c). 
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Figure 13A contains colored circles representing the re-calculations intended to match the original curve 
A. Each circle was obtained by selecting a value of r a and executing the NACA 496 solution method for 
Subcase 1. The agreement between the original and present results is very good at values of r a of 0.75 
and greater, but slightly less so at lower values of r a . If original curve A is extrapolated to the left, the 
difference between the original and present values of Fatr a = 0.25 is about 7.5%. The large X's at 
values of r a of 0.50, 0.75, and 1.00 are from the independent calculation and are in excellent agreement 
with the present calculations. 

Figure 13B(a) contains colored squares representing the re-calculations intended to match the original 
curve B. Each square was also obtained by selecting a value of r a and executing the NACA 496 solution 
method for Subcase 1. The trends of the original and present results agree (F decreases with increasing 
values of r a ), but in terms of absolute numbers the agreement is very poor: at the higher values of r a , 
the magnitudes of the original results are more than a factor of two larger than those of the present 
results; at the lower values of /^they approach a factor of two smaller. In fact, the vertical scale of the 
original figure had to be extended in order to accommodate the present result at r a = 0.25. 

N.B.: In the following discussion there are two uses of the word "square." The first use refers to the 
product of a quantity and itself; the second use refers to the geometric shape of the colored plotting 
symbols in figure 13B. 

In contemplating the poor agreement between the original and present results, the present author 
noticed that the values of the present results, represented by the colored squares, at r a = 0.375 and 1.0 
(about 1.5 and 0.25, respectively) are approximately equal to the squares of the values of the original 
curve at the same r a ’s (about 1.25 and 0.5, respectively). This relationship led the present author to 
compute and plot the squares of the values of the original curve B (curve labeled in the figure as the 
"square of curve B"). 

As can be seen in figure 13B(b), the agreement between the present results and the square of curve B is 
much better than the agreement between the present results and curve B itself. The better agreement 
may be quantified by examining the sums of the squared differences between the results: those of the 
former are less than half those of the latter. 

This better agreement suggests that an error occurred in the creation of the original curve B. Because 
the author of NACA 496 uses the symbol Fto mean, at times, a quantity and, at other times, the square 
of that quantity, there is ample room for one meaning to be used when the other meaning was 
intended. The present author believes such an error occurred in the creation of curve B (but not in the 
creation of curve A). 

In figure 13B(c), the large X's at values of r a of 0.50, 0.75, and 1.00 are from the independent calculation 
and are in excellent agreement with the present calculations and is supporting evidence that the original 
curve B was incorrectly plotted. 
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Variation of Center of Gravity Location . - The original figure 14 illustrates the dependency of F on the 
location of the center of gravity, x a , for three different combinations of quantities. Each combination is 
represented by a different curve, labeled A, B, and C. According to the sign convention established in 
figure 2 of N AC A 496, increasingly positive values of x a corresponds to an aft movement of the center of 
gravity. 

Because of the complexity of the results to be presented for this figure and in an attempt to preserve 
the figure numbering convention described at the beginning of this section, figure 14 will be divided up 
into "offspring" figures: 14A, 14B, and 14C. In addition, each offspring figure will have parts (a), (b), and 
(c). 

According to the figure legend in the original figure 14, the quantities specific to this example are 

for curve A: k-V^ oo; a = -0.4; coi/(o 2 - Vs) r a - 0.5; x a variable 

and 

for curve B: k=V i; a = -0.4; coi/co 2 - 34; r a - 0.5; x a variable 

and 

for curve C: k-Y 1 00 ; a = -0.4; (Oi/(o 2 - 1; r a = 0.5; x a variable. 

The colored circles in Figures 14A are the result of present calculations intended to match the original 
curve A; the colored squares in Figures 14B, to match the original curve B; and the colored diamonds in 
Figures 14C, to match the original curve C. Each present result represented by a symbol was obtained 
by selecting a value of x« and executing the NACA 496 solution method for Subcase 1. 

The original figure 14 contains conflicting information regarding curve C. The figure legend defines the 
quantities listed above; but, within the figure itself is a notation indicating other quantities - "coi/cd 2 = 1, 
4 k- .01" - with a leader pointing to curve C. The ratio of frequencies is the same as the ratio specified 
in the figure legend; however, the value of >400) differs from that in the figure legend. Present 
calculations were performed for both values of k. 

In figure 14A(a), the original results and the present results differ significantly, showing opposite trends 
with increasing values of x a - However, the present author observed the following: the present value of 
F at x a - 0 is very close to the original value of F at x a = 0.4 and vice versa; the present value of F at x a = 
0.1 is very close to the original value of F at x a = 0.3 and vice versa; and the present value of F at x a - 0.2 
is very close to the original value of F at x a = 0.2. 

This observation led to figure 14A(b), in which the original curve A has been rotated about the vertical 
axis x a = 0.2. The rotated curve A and the colored circles exhibit excellent agreement. 

In Figure 14A(c) the large X's at values of x a of 0.1, 0.2, and 0.3 are from the independent calculation and 
are in excellent agreement with the present results and the rotated curve A. This triple agreement 
suggested to the present author that original curve A was incorrectly plotted. 
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In the early 1930's, intricate and multi-stepped engineering calculations (such as those required to 
produce the results in NACA 496) were performed by hand; results were recorded and transcribed by 
hand; and then plotted by hand. It is not unrealistic to consider that with the many steps performed by 
humans, errors occurred. The present author believes that in the case of curve A, transcription errors 
occurred - correct values of F were paired with non-corresponding values of x a - which then led to an 
incorrect plot. 

The trends observed in figures 14A are repeated in figures 14B and 14C: in figures 14B(a) and 14C(a) the 
original results and the present results (represented now by the colored squares and diamonds) differ 
significantly, showing opposite trends with increasing values of x a ; in figures 14B(b) and 14C(b) the 
original curves B and C have been rotated about the vertical axis x a = 0.2 producing excellent agreement 
with the present results; and in figures 14B(c) and 14C(c) the large X's from the independent calculation 
are in excellent agreement with the present results and the rotated curves B and C. Again, this time for 
curves B and C, the present author believes that transcription errors occurred - correct values of F were 
paired with non-corresponding values of x a - which then led to incorrect plots. 


Wings A and B 

The final three figures in NACA 496 (figs. 15, 16, and 17) present results for two wings, designated in 
NACA 496 as Wing A and Wing B. For each wing, experimental results were obtained. Both wings had 
symmetrical airfoils with chords of 12.7 centimeters (5 inches) and thickness-to-chord ratios of 0.12. 
Both wings were tested at zero degrees incidence. 

The text of Appendix II of NACA 496 states that Wing A was constructed of aluminum with quantities: 

K— V 416 ; a = -0.4; x a = 0.31, 0.173, and 0.038; r a 2 - 0.33; C0a=7x2n. 

The text of Appendix II of NACA 496 states that Wing B was constructed of wood, with a trailing-edge 
control surface, with quantities: 

fc=Y\oo] c = 0.5; a = -0.4; x a - 0.192; r a 2 - 0.178; 

Xjg= 0.019; r/ = 0.0079; <%= 17.6x271. 

Figure 15 contains Subcase-1 results for Wing A in which flutter velocity is plotted as a function of 
frequency ratio coh/coa . 

Figure 16 contains Subcase-2 results for Wing B in which flutter velocity is plotted as a function of 
frequency ratio cop/con . 

Figure 17 contains Subcase-3 results for Wing B in which flutter velocity is plotted as a function of 
frequency ratio co a /cop. 
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To obtain the curves in these figures, first the appropriate subcase-specific NACA 496 solution method 
was executed, resulting in curves of subcase-specific F (as defined) vs. subcase-specific Q. Next, for each 
value of F equation (13) was employed, converting the ordinates in these curves from nondimensional F 
to dimensional v/. Finally, the equations for the subcase-specific £)’s from equations (5a), (5b), and (5c) 
were employed, converting the abscissas from subcase-specific Q to subcase-specific frequency ratio. 

Referring back to equation (13), it can be seen that to obtain the dimensional flutter velocity reference 
values co r and r r are required, which differ from subcase to subcase. For Subcase 1, values of co a and r a 
are required; for Subcase 2, a value of a>h is required; and for Subcase 3, values of copsnd rp are 
required. NACA 496 supplies values for co a , r a , and rp, but, in another tripping point, it does not supply 
values for coh and cop. Therefore, for the present calculations for Subcases 2 and 3, values of coh and cop 
had to be chosen. The choice was based on a trial-and-error investigation of candidate values and the 
selection criterion for the chosen values of (Dh and cop was that there be reasonable agreement between 
the original and the present curves. 

Subcase 1 . - Figure 15 contains theoretical and experimental results, in black. The present results are 
shown in colored O and + symbols and solid and dashed colored curves. Even though three values of x a 
are listed above in the quantities for this subcase, x a = 0.173 was chosen for the present calculation 
because it produced the best agreement with the original curve in figure 15. Over the entire frequency 
range the trends for original theoretical results and both sets of present results are in very good 
agreement. The overall agreement over the first half of the range is excellent with maximum 
differences less than 5 percent; over the second half, less so with maximum differences of about 15 
percent. 

Subcase 2 . - Figure 16 contains theoretical and experimental results, in black. The present results are 
shown in colored O and + symbols and solid and dashed colored curves. Theoretical results are given 
for two values of quantity x^: xp- 0.019 and xp- 0.01. For this subcase, an assumed value of con of 
5.8x2tt was selected because it produced reasonable agreement between the original results and the 
present results. 

In examining the colored curves, it is seen that the upper (corresponding to higher values of v) and 
lower (lower values of v) extremes of both curves extend to a frequency ratio of zero. Both colored 
curves are comprised of solid and dashed portions, indicating that two roots of equation (10) are 
present in each solution. The transition from solid to dashed occurs between frequency ratios of 0.8 
and 0.9. The calculations that produced these curves had a range of 1/k from 0 to 500. 

In examining the black theoretical curves, it is seen that the upper extremes of both curves terminate at 
a frequency ratio of about 0.75 and the lower extremes terminate at a frequency ratio of about 0.82. 

The difference in horizontal extent between the black and colored curves suggests to the present author 
that the calculations that produced the original curves had a very limited range of 1/k. Further, because 
the frequency ratio at the lower extremes of the black curves corresponds roughly to the frequency 
ratio at the transition of the colored curves from solid to dashed, it is possible that the black curves only 
capture a single root of equation (13/XXVI). 
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Subcase 3 . - Figure 17 contains theoretical and experimental results, in black. The present results are 
shown in colored O and + symbols and solid and dashed colored curves. For this subcase, an assumed 
value of (Qpo\ 4.4x271 was selected because it produced reasonable agreement between the original 
results and the present results. 

As was observed in figure 16 for subcase 2 because the frequency ratio at the lower extreme of the 
black curve corresponds roughly to the frequency ratio at the transition of the colored curve from solid 
to dashed, it is possible that the black curve only captures a single root of equation (13/XXV). 


VII. CONCLUDING REMARKS 

This paper has reported on an extensive study of reference 2, NACA Report No. 496, "General Theory of 
Aerodynamic Instability and the Mechanism of Flutter," by Theodore Theodorsen (referred to herein as 
NACA 496), including the re-computation of its numerical results. The study included reading and re- 
reading the paper, checking its equations, and checking its solution methods. With the few exceptions 
noted in the present paper, the study revealed the equations and solution methods of NACA 496 to be 
correct. Because of the crude (by today's standards) computational tools available at the time NACA 
496 was written, the author of NACA 496 conceived of and employed several time- and effort-saving 
computational shortcuts, pointed out and summarized in a table in the present paper. The study also 
revealed some "tripping points," aspects of NACA 496 that have the potential to cause confusion, and 
identified some errors in NACA 496 that were not corrected in the later versions of the paper. Tripping 
points and errors are also summarized in tables in the present paper. 

Appendix II in NACA 496 contains a number of numerical examples, all of which were re-computed in 
the present paper. The re-computed results were overlaid on the original figures from NACA 496. The 
re-computations were accomplished by employing the method developed in the original work, but by 
using modern computational tools. In most instances the overall shapes of the original and re- 
computed curves are similar, with varying degrees of agreement in magnitude. In some instances, there 
is significant dissimilarity between the original and re-computed curves. Most dissimilarities are 
attributed to human computational and transcription errors that are speculated to have occurred during 
the creation of NACA 496. One dissimilarity (namely, that regarding the quantity F in fig. 12) is 
unattributed. Several of the re-computed results were checked by independent calculations and, in 
each instance, the re-computed and independent results were in excellent agreement. 

Appendix A of the present paper contains NACA 496. 

Appendix B contains a Matlab® program that implements the two-degree-of-freedom flutter equations 
and the two-degree-of-freedom solution method found in Appendix I of NACA 496. By executing this 
program, results in figures 5 through 11 and 15 through 17, found in Appendix II of NACA 496, may be 
re-computed directly. This program produces both physical and non-physical solutions. Physical 
solutions are based on real values of the quantity Fand the /7s; non-physical solutions are based on 
complex values. 
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Appendix C of the present paper offers three alternate solution methods, all of which derive from the 
two-degree-of-freedom flutter equations of NACA 496. Compared to the number of computations 
involved in the NACA 496 solution method, these alternate solution methods are much more 
computationally intensive and would have been excluded from consideration in 1934 for that reason. 

To accuracies within a tenth of a percent, all three alternate solution methods yielded the same value of 
flutter velocity, v/, and the same value of flutter reduced frequency, kf, as those predicted by the 
solution method of NACA 496. 

Appendix D of the present paper offers a solution of the three-degree-of-freedom flutter equations 
following the method proposed in NACA 496. This solution method produced an example result for the 
"standard case." In addition, for systematic variations in the dimensional modal frequencies, in the 
limit, the three-degree-of-freedom results (v/and kf) approached the two-degree-of-freedom results, 
including the number of flutter modes predicted. 
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AERONAUTIC SYMBOLS 

1. FUNDAMENTAL AND DERIVED UNITS 


w 

g 

m 

I 

P 

S 

S„ 

G 

b 

c 

A 

V 

<z 

L 

D 

D 0 

D t 

D v 

C 



Symbol 

Metric 

English 

Unit 

Abbrevia- 

tion 

Unit 

Abbrevia- 

tion 

Length 

l 

meter 

m 

foot (or mile) _ 

ft (or mi) 

Time 

t 

second _ _ 

s 

second (or hour) 

sec (or hr) 

Force. __ . 

F 

weight of 1 kilogram 

kg 

weight of 1 pound-- 

lb 

Power 

P 

horsepower (metric) 


horsepower 

hp 

Speed 


/kilometers per hour 

kph 

miles per hour 

mph 

V 

(meters per second 

mps 

feet per second 

fps 


2. GENERAL SYMBOLS 


Weight =mg 

Standard acceleration of gravity =9.80665 m/s 2 
or 32.1740 ft/sec 2 


Mass=— 

9 

Moment of inertia =ra& 2 . (Indicate axis of 
radius of gyration k by proper subscript.) 
Coefficient of viscosity 


v Kinematic viscosity 

p Density (mass per unit volume) 

Standard density of dry air, 0.12497 kg-m -4 -s 2 at 15° C 
and 760 mm; or 0.002378 lb-ft -4 sec 2 
Specific weight of “standard” air, 1.2255 kg/m 3 or 
0.07651 lb/cu ft 


3. AERODYNAMIC SYMBOLS 


Area 

Area of wing 
Gap 
Span 
Chord 

Aspect ratio, ^ 

True air speed 
Dynamic pressure, ipF 2 


Lift, absolute coefficient C L = 


Drag, absolute coefficient C D = 


qS 

D 


qS 


D„ 


Profile drag, absolute coefficient D Dq ~-^ 


D t 


Induced drag, absolute coefficient 


D„ 


Parasite drag, absolute coefficient 0^=^ 

C 

Cross-wind force, absolute coefficient 
2626° 


i w Angle of setting of wings (relative to thrust line) 
i t Angle of stabilizer setting (relative to thrust 
line) 

Q Resultant moment 

ft Resultant angular velocity 

VI 

R Reynolds number, p — where l is a linear dimen- 
sion (e.g., for an airfoil of 1.0 ft chord, 100 mph, 
standard pressure at 15° C, the corresponding 
Reynolds number is 935,400; or for an airfoil 
of 1.0 m chord, 100 mps, the corresponding 
Reynolds number is 6,865,000) 
a Angle of attack 

€ Angle of downwash 

a 0 Angle of attack, infinite aspect ratio 

at Angle of attack, induced 

a a Angle of attack, absolute (measured from zero- 

lift position) 

7 Flight-path angle 
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GENERAL THEORY OF AERODYNAMIC INSTABILITY AND THE MECHANISM OF 

FLUTTER 


By Theodore Theodorsen 


SUMMARY 

The aerodynamic jorces on an oscillating airfoil or 
airfoil-aileron combination oj thxee independent degrees 
of freedom hare been determined. The problem resolves 
itself into the solution of certain definite integrals , which 
have been identified as Bessel functions of the first and 
second kind and of zero and first order. The theory, 
being based on potential flow and the Kutta condition, 
is fundamentally equivalent to the conventional wing- 
section theory relating to the steady case. 

The air forces being known , the mechanism of aerody- 
namic instability has been analyzed in detail. An exact 
solution, involving potential flow and the adoption of the 
Kutta condition, has been arrived at. The solution is of 
a simple form and is expressed by means of an auxiliary 
parameter k. The mathematical treatment also provides 
a convenient cyclic arrangement permitting a uniform 
treatment of all subcases of two degrees of freedom. The 
flutter velocity, defined as the air velocity at which flutter j 
starts, and which is treated as the unknown quantity, is j 
determined as a function of a certain ratio of the fre- ! 
quencies in the separate degrees of freedom for any magni- 
tudes and combinations of the airfoil-aileron parameters. 

For those interested solely or particularly in the numeri- 
cal solutions Appendix I has been prepared. The rou- 
tine procedure in solving numerical examples is put 
down detached from the theoretical background of the 
paper. It first is necessary to determine a certain number 
of constants pertaining to the case, then to perform a few 
routine calculations as indicated. The result is readily 
obtained in the form of a plot of flutter velocity against 
frequency for any values of the other parameters chosen. 
The numerical work of calculating the. constants is sim- 
plified by referring to a number of tables , which are in- 
cluded in Appendix I. A number of illustrative examples 
and experimental results are given in Appendix II. 

INTRODUCTION 

It has been known that a wing or wing-aileron struc- 
turally restrained to a certain position of equilibrium 
becomes unstable under certain conditions. At least 
two degrees of freedom are required to create a con- 
dition of instability, as it can be shown that vibrations 


of a single degree of freedom would be damped out by 
the air forces. The air forces, defined as the forces due 
to the air pressure acting on the wing or wing-aileron 
in an arbitrary oscillatory motion of several degrees of 
freedom, are in this paper treated on the basis of the 
theoiy of nonstationary potential flow. A wing- 
section theory and, by analogy, a wing theory shall be 
thus developed that applies to the case of oscillatory 
motion, not only of the wing as a whole but also to 
that of an aileron. It is of considerable importance 
that large oscillations may be neglected ; in fact, only 
infinitely small oscillations about the position of 
equilibrium need be considered. Large oscillations 
are of no interest since the sole attempt is to specify 
one or more conditions of instability. Indeed, no 
particular type or shape of airfoil shall be of concern, 
the treatment being restricted to primary effects. The 
differential equations for the several degrees of freedom 
will be put down. Each of these equations contains a 
statement regarding the equilibrium of a system of 
forces. The forces are of three kinds: (1) The inertia 
forces, (2) the restraining forces, and (3) the air forces. 

There is presumably no necessity of solving a general 
case of damped or divergent motion, but only the 
border case of a pure sinusoidal motion, applying to the 
case of unstable equilibrium. This restriction is par- 
ticularly important as the expressions for the air force 
developed for oscillatory motion can thus be employed. 
Imagine a case that is unstable to a very slight degree; 
the amplitudes will then increase very slowly and the 
expressions developed for the air forces will be appli- 
cable. It is of interest simply to know under what 
circumstances this condition may obtain and cases in 
which the amplitudes are decreasing or increasing at a 
finite rate need not be treated or specified. Although 
it is possible to treat the latter cases, they are of no 
concern in the present problem. Nor is the internal 
or solid friction of the structure of primary concern. 
The fortunate situation exists that the effect of the 
solid friction is favorable. Knowledge is desired con- 
cerning the condition as existing in the absence of the 
internal friction, as this case constitutes a sort of lower 
limit, which it is not always desirable to exceed. 
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Owing to the rather extensive field covered in the 
paper it has been considered necessary to omit many 
elementary proofs, it being left to the reader to verify 
certain specific statements. In the first part of the 
paper, the velocity potentials due to the flow around 
the airfoil-aileron are developed. These potentials 
are treated in two classes: The noncirculating flow 
potentials, and those due to the surface of discon- 
tinuity behind the wing, referred to as “circulatory ” 
potentials. The magnitude of the circulation for an 
oscillating wing-aileron is determined next. The 



Figure 1.— Conformal representation of the wing profile by a circle. 

forces and moments acting on the airfoil are then 
obtained by integration. In the latter part of the 
paper the differential equations of motion are put 
down and the particular and important case of un- 
stable equilibrium is treated in detail. The solution 
of the problem of determining the flutter speed is 
finally given in the form of an equation expressing a 
relationship between the various parameters. The 
three subcases of two degrees of freedom are treated 
in detail. 

The paper proposes to disclose the basic nature of 
the mechanism of flutter, leaving modifications of the 
primary results by secondary effects for future investi- 
gations. 1 Such secondary effects are: The effects of a 
finite span, of section shape, of deviations from poten- 
tial flow, including also modifications of results to 
include twisting and bending of actual wing sections 
instead of pure torsion and deflection as considered in 
this paper. 

The supplementary experimental work included in 
Appendix II similarly refers to well-defined elementary 
cases, the w r ing employed being of a large aspect ratio, 
nondeformable, and given definite degrees of freedom 
by a supporting mechanism, with external springs 
maintaining the equilibrium positions of wing or wdng- 
aileron. The experimental w r ork w r as carried on 
largely to verify the general shape of and the approxi- 
mate magnitudes involved in the theoretically pre- 
dicted response characteristics. As the present report 
is limited to the mathematical aspects of the flutter 
problem, specific recommendations in regard to prac- 
tical applications are not given in this paper. 

1 The effect of internal friction is in some cases essential; this subject will be 

contained in a subsequent paper. 


VELOCITY POTENTIALS, FORCES, AND MOMENTS OF 
THE NONCIRCULATORY FLOW 

We shall proceed to calculate the various velocity 
potentials due to position and velocity of the individ- 
ual parts in the whole of the wing-aileron system. 
Let us temporarily represent the wing by a circle (fig. 
1). The potential of a source c at the origin is given 
by 

v?=^l°g (x 2 + y 2 ) 

For a source c at {x u y x ) on the circle 

(p= tn - log { (z — #i) 2 + (y~Vi) 2 } 


Putting a double source 2c at {x x ,y x ) and a double 
negative source —2c at (x lf — y x ) we obtain for the flow 
around the circle 


*=^i° g 


(x-x^+jy-yj 2 

(x-x^+iy + yi ) 2 


The function <p on the circle gives directly the sur- 
face potential of a straight line pq, the projection of the 
circle on the horizontal diameter. (See fig. 1.) In 
this case y = -y / 1 — x 2 and ^ is a function of x only. 

We shall need the integrals: 


s. 


and 




I 


log 


(x—x\) 2 + (y—yi) 2 

(x—xi) 2 + (y+2/1) 2 


fo —c)<\x x = — -yjl—C 2 -y]l —X 2 


where 


- cos -1 c(z — 2c) Vl — x 2 + (x — c) 2 \ogN 
N 1 — cx— Vl — a^Vl — c 2 


The location of the center of gravity of the wing- 
aileron x a is measured from a , the coordinate of the 
axis of rotation (fig. 2); xp the location of the center 



Figure 2.— Parameters of the airfoil-aileron combination. 


of gravity of the aileron is measured from c, the coordi- 
nate of the hinge; and r a and rp are the radii of gyration 
of the wing-aileron referred to a , and of the aileron 
referred to the hinge. The quantities xp and rp are 
“reduced” values, as defined later in the paper. The 
quantities a, x a , c, and Xp are positive toward the rear 
(right), h is the vertical coordinate of the axis of rota- 
tion at a with respect to a fixed reference frame and is 
positive downward. The angles a and /3 are positive 
clockwdse (right-hand turn). The wind velocity v is to 
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the right and horizontal. The angle (of attack) ot 
refers to the direction of v, the aileron angle p refers to 
the undeflected position and not to the wind direction. 
The quantities r a and r 0 always occur as squares. 
Observe that the leading edge is located at —1, the 
trailing edge at +1. The quantities a, c } x a , xp, r a , 
and which are repeatedly used in the following 
treatment, are all dimensionless with the half chord b 
as reference unit. 

The effect of a flap bent down at an angle p (see fig. 
2) is seen to give rise to a function ip obtained by sub- 
stituting — vpb for e; hence 

<p fi = [*s/l—x 2 cos-’c -- (x - c) log A r ] 

To obtain the effect of the flap going down at an 
angular velocity p , we put * = - fa — c)pb 2 and get 

n - g-[ Vl - c 2 Vl - * + cos“‘c(z- 2c) 

— (x — c) 2 log N] 

To obtain the effect of an angle a of the entire air- 
foil, we put c= — 1 in the expression for <p p , hence 
(p a ~vab^l — x 2 

To depict the airfoil in downward motion with a veloc- 
ity h (+ down), we need only introduce ^ instead of a. 

Thus . 

<p- h = kb^l—x 2 

Finally, to describe a rotation around point a at an 
angular velocity a, we notice that this motion may be 
taken to consist of a rotation around the leading edge 
c — — 1 at an angular velocity a plus a vertical motion 
with a velocity — a (1 + a) b. Then 

<pi = ^-ir(x + 2) Vl —x 2 —a(l -f a)b 2 -yjl — x 2 

2ir 

— a& 2 ^~ x — Vi “ 

The following tables give in succession the velocity 
potentials and a set of integrals 2 with associated con- 
stants, w r hich we wdll need in the calculation of the air 
forces and moments. 

VELOCITY POTENTIALS 

(Pa^vaby 1 1 —x 2 
(pit — hb \ 1 x 2 

1— IE 2 cos ~‘c—(x — c) log AH 

7T 

^8 = ^- $b 2 [\' 1 — C 2 -y/ 1 — X 2 + (x — 2 c)^\—X 2 COS'^ 

2i r 

— (x — c) 2 log N] 

where 1 -cx~ Vl -^Vl ~ c 2 

x — c 

1 Some of the more difficult integral evaluations are given in Appendix III. 


<p a dx — — -yVa.T± 

^dx = — 2^ 4 
<Pa dx = ab 2 l\ 


vpT, 
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b 2 • 

iP8^X— — 


j: 
j: 
r 

J <Ppdx= - 

i 

tp a (x— c)dx= — yfiaTi 
J' (Ph(x-c)dx= - ^hTi 
(pi(x — c)dx = ab 2 T ld 
f <p8(x — c)dx= —~vpT 2 

Jc xr 

X i b 2 . 

ipfa-c) dx= ~ 2 -PT* 


INTEGRALS 
'•+1 


<Padx = 2 Vair 

+1 , b s 

<Phdx — o n,7r 


,dx = — ab 2 


r: 

i: 
j> 

r«d.— i 

I> 

/> 

r+i 

I <ph(x-c)dx = 

s: 


>vPT 4 • 


A 2 . 

idx = - 2 ^ 


+i b 

ip a (x c) dx = — ^vacir 


- ^hcir 


<Pa (x ~c)dx= ab 2 T u ir 

b „ 
PP % 8 


r + 1 j 

I <pfi(x-c)dx= -- 

r+i 

I ^ <p(){x — c)dx = 


b 2 pT 7 

2 


CON.STANTS 


Ti = — g Vi “ c 2 (2 + c 2 ) + c cos _1 c 

T 2 = c(1 — c 2 ) — Vl — c 2 (l -fc 2 )cos -1 c + c(cos -1 c) 2 

T z = - (| + c 2 ) (cos' ‘c) 2 + \c -v/l 3 ? cos-'c(7 + 2c 2 ) 


8 


(1 — c 2 ) (5c 2 +4) 


1\ = — COS _1 C + cVl —C 2 

T5 — ■ (1 —c 2 ) — (cos^c) 2 \- 2c ijl — c 2 cos ~ l c 

T t =T 2 

T, = - (| + c 2 ) cos-'c + g c VI 3 ? (7 + 2c 2 ) 

T s — — 5 Vl — ? (2c 2 + 1) + c cos -1 c 
o 

r, = | [| (■ Vl 3 ?) + aT 4 ] = \ (-P + aT t ) 

where p — — |(vi — c 2 ) 3 
T, 0 — Vl — c 2 + cos' 1 c 
T n = cos" 1 c (1 - 2c) + Vl 3 ? (2 - e) 
ru-vr 3 ? (2 + c) — cos -1 c (2c + 1 ) 

y T i3 = |[-^ 7-(c-a)T 1 l 

rp 1,1 

T h 16 + 2 ac 

FORCES AND MOMENTS 

The velocity potentials being known, we are able to 
calculate local pressures and by integration to obtain 
the forces and moments acting on the airfoil and 
aileron. 
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Employing the extended Bernoulli Theorem for un- 
steady flow, the local pressure is, except for a constant 



where w is the local velocity and p the velocity poten- 
tial at the point. Substituting w= ^ v + ^j~ we obtain 

ultimately for the pressure difference between the 
upper and lower surface at £ ^ 

where v is the constant velocity of the fluid relative to 
the airfoil at infinity. Putting down the integrals for 
the force on the entire airfoil, the moment on the flap 


y 



Figure 3.— Conformal representation of the wing profile with reference to the 
circulatory flow. 


around the hinge, and the moment on the entire air- 
foil, w r e obtain by means of partial integrations 
r + 1 

P=—2pb\ v>d^ 

Mp = — 2pb 2 £ <p(x — c) dx + 2pvbj <pdx 

M a = — 2pb 2 j* ^ (p(x — c)dx + 2pvb^ <pdx 
r+i 

— 2pb 2 I ^ 4> (c - a) dx 

Or, on introducing the individual vclocit} 7 potentials 
from page 5, 

P= — pb 2 [vtt6l-\- irh — bwaa— vT 4 (3 — bT$] ( 1 ) 

Mp = — pb* — vT } a — 7\h + 2 I\ 3 bd — ^ vT 2 (3 — ^ 

+ pvb 2 jj - v'T.a - T t h + 2TJ)a - ^ vTrf - i 7’ 2 6/3 J 
~ — pb~ I” 7\v~a — (27g 4 T,) boa 4 2 r l\-J) 2 a 4 — T 5 v 2 ($ 

L 7T 

+ T bv& - ± b 2 T,& + T<vK - 7 ’,&/;] (II) 

M a — — pb 2 £ — TV 2 a + 7 + a 2 ^ b 2 oc + v 2 T 4 p + { T x - T 
- (c - a) T 4 } bfiv+ {- T 7 -(c- a) T x }b 2 $ 
bunk— (III) 


VELOCITY POTENTIALS, FORCES, AND MOMENTS 
OF THE CIRCULATORY FLOW 


In the following we shall determine the velocity 
potentials and associated forces and moments due to a 
surface of discontinuity of strength TJ extending along 
the positive x axis from the wing to infinity. The 
velocity potential of the flow around the circle (fig. 3) 
resulting from the vortex element - Ar at (X 0 , 0) is 


* ■ t? ~ tan_1 ^i i j 


AT 

= 2tt 


tan" 




X 2 


0 Yo+ i„) 


^ *4 l' 2 4-1 


where ( X , Y) are the coordinates of the variable 
and X 0 is the coordinate of — AF on the x axis. 


Introducing A' 0 + -^ = 2 j ( , 

'*0 

or X 0 = x 0 4- V^o 2 “ 1 on the x axis 
and X~x and T=yl — .r 2 on the circle 
the equation becomes 


A1 . 

^zx«= - tail' 


-i Vi -r y ¥ 


1 — xx Q 

This expression gives the clockwise circulation 
around the airfoil due to the element -AF at x 0 . 


We have:p=-2p(^4»g) 


But, since the element — AF will now be regarded as 
moving to the right relative to the airfoil with a 
velocity v 

dip __ d<p 

dt ~ dx 0 v 



F urther 

* , ^ o yi - x 2 

2n djp _ / , (1 -xx 0 ) Vl -y 2 (l-yj 0 ) 2 

AF d* v<r ° , , (l-rOW-l) 

(l-^o) 2 

VVEI 1 

Vl - x 2 {jCq-x) 

and 

1 Xo /~2 ~ 7 JL 

2*$*^ r — -^ {\—xxq) \Xq—\ 1 X ° (1-xxp) 2 

Al'd-fo V . , (l-jflfcr,*-!) 

(1 - xx 0 ) 2 

V 1 — X 2 1 

W I (*o - x) 

By addition: 

dip dip = AT x 0 Yx 

dx dx 0 ~ 2 w Vr^VV 3 ^ 
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Finally 


To obtain the force on the aileron, we need the 
integral 

Xp + X 




■yjx o 2 — 1 V 1 


dx 


ATf *o _i , V I _ z 2 "T 

2-V^l C0S X + lW=li 

, 2 -iJ 


= ATf Xp 


, VI -e 

= COS 'c + J 7 = 

1 Vxo‘ 


'2xLVi7 

Thus, for the force on the aileron 


AP, 


= cos -, c H — j- 
1 V^o‘ 




os"‘c— VW) 


- 

+ VP^] 

Integrated, with Ar= Udx 0 
P ‘* = -^[(cos-c- 


+ Vi- 




for c= — 1 we obtain the expression for P, the force 
on the whole airfoil 




(IV) 


Since U is considered stationary with respect to the 
fluid elements 

U=i(vt-x 0 ) 

where t is the time since the beginning of the motion. 
U is thus a function of the distance from the location 
of the first vortex element or, referred to a system 
moving with the fluid, U is stationary in value. 
Similarly we obtain for the moment on the aileron 

P/d^ j^V AT P (z-e) (*, + ») 

JAe>* ’ 2xJ c vr^v'xT^r 


dz 


= '"^^Tt oVl “ l2+?3 ^" cVl “ a:2 

+ Q-ZoC^cos _1 :eJ 

= + fiv^[ ( * 0+ ^ c)VI ^ 

+ - 2 x 0 c^cos~ , cJ 

= + 2VLVi7^ ( - vl_c ccos c > 

• 1 

(cos'" *c - c Vl - C 2 )J 




4M ,__„^[^|vr=e( 1+ «) 

Putting AT = L 7 dr 0 and integrating 


— cos -1 c^ 


, c+; 


'i V^o 2 -l 


- h dx 0 


+ (c 


COS L C — cV (1 — c* 

Further, for the moment on the entire airfoil around a 

S-Xt + t) (x ~ a > Ax = - tv vii[(*“ + i“ a ) vrr? 

+ G-^) cos " i ]-! =+ t!r; 7 xViG” x » o > 


Hf VSH (V) 


and AM a — — pvb 2 A V 

Integrated, this becomes 

1 

2 x <& 

V^T^ 0 


1 

2 ^ 0 ® 

JxJ~-\ 


M a = — 

= -pvb 2 
= - pvb 2 


I 


| 1,1 

2 + 2 X ° X{ 


<" 4 D 


I/djo 


V -Co 2 1 V-Co 2 1 

1 <»■> 


X n - 


THE MAGNITUDE OF THE CIRCULATION 

The magnitude of the circulation is determined by 
the Kutta condition, which requires that no infinite 
velocities exist at the trailing edge, 
or, at x = 1 

d - tpa F (ph "t" <pa F (pfi 4" ipp) = finite 
Introducing the values of <p a , etc. from page 5 and 


, from 


^ page 6 gives the important r< 


•elation: 


i r-v 

2ttJi y 


V®o± 1 

Jx 0 — 1 


Udx 0 = va+h + b\ 


h+bQ- a y 


+‘^»p+bp i fi 

ir w 7 r 


(VII) 


This relation must be satisfied to comply with the 
Kutta condition, which states that the flow shall leave 
the airfoil at the trailing edge. 

It is observed that the relation reduces to that of the 
Kutta condition for stationary flow on putting ar 0 = 00 , 
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and in subsequence omitting the variable parameters 
a, $, and 
Let us write 


Introduced in (IV) 

P— —2t pvbQ 

from (V) 

Mfi — — 2pvb 2 

j; 


f°° *0 

Jl Vzo 2 — 1 


Udx 0 


(vi =?(l + §) - cos- C (c + 5)): 


X 


Xq 


V^O 2 - 1 


Udx 0 


nm 


Udx ( 


+ ^COS -1 C — cVl “ C 2 ^ 



M« = — 2irpvb 2 
Introducing 


p a+ i) 


r_ xp 

l\Ji Vzp 2 - 1 


Udx, 


r l x 0+1 

J Vzo-l 


£/dz 0 


<7= 




rC/da-o 


we obtain finally 


f" lx o+l 

Ji ViFT 

P=-2pvbirCQ (VIII) 

A/j = — 2p»4 2 [(vi -c 2 (l + |) - cos -1 c(c + |))(7 

+ § (cos- c-cV1-c 2 )]q=-p»6 2 (7’i 2 C- r 4 )Q (IX) 
= 2 Trprft 2 [(a + D^-ljQ (X) 


where Q is given above and C= C(k ) will be treated in 
the following section. 

VALUE OF THE FUNCTION C(Jt) 

put ?j= tv ‘[‘(KM 

where s = vt (s— >a >), the distance from the Jirst vortex 
element to the airfoil, and k a positive constant deter- 
mining the wave length, 
then 


C(k) = 


r 


Xq 


s: 


V*o 2 - 1 


e ikx odx. 


x ° + ie-«'odx„ 


(XI) 


V^T 

These integrals are known, see next part, formulas 
(XIV)— (XVII) and we obtain 3 


C(k)=- 


i 


-J t + i¥t 


7T J _^y I’][y _A?L T (^ld" F<])+i(Fi Jo) 
2^1 

( ■ - Jl +iF,)[ - (J, + Yo) - i(F , - Jo)l 

(j,+y„) 2 +(F,-j„) 2 

J,(J,+ F 0 )+F 1 (F 1 -J 0 ) 

(J 1 +F„) 2 +(F,-J „) 2 

. F 1 ( t /,+ F 0 )-J,(F,-Jo) _ 

1 (J 1 +Fo) 2 +(F 1 -J 0 ) s 

where 

p e/ 1 (J 1 +Fo)+F 1 (F 1 -Jo) 

(J r i+F 0 ) 2 +(F 1 -e/ 0 ) 2 CAA1) 


G= — 7 


FiFq + Jic/q 


(J 1 +F„) 2 +(F 1 - J „) 2 (XIII) 

These functions, which are of fundamental import- 
ance in the theory of the oscillating airfoil are given 

graphically against the argument ^ in figure 4. 

SOLUTION OF THE DEFINITE INTEGRALS IN C BY MEANS OF BESSEL 
FUNCTIONS 

We have 


K r 


-«-f 


g-*COBht cog k ^ 


where 


and 


but 


(Formula (34), p. 51— Gray, Mathews 
& MaeRobert: Treatise on Bessel 
Functions. London, 1922) 

in v 

K, (<) = e 2 GJit) 

(Eq. (28), sec. 3, p. 23, same reference) 
Q. W ! - - Yn M + [log 2 - y + (X) 

F»C*0=f Y n (x) + (log 2-y)J n (x) 

(where Y n (x) is from N. Nielsen: 
Handbuch der Theorie der Cylinder- 
funktionen. Leipzig, 1904). 


* This may also be expressed in Hankel functions, C— 

(V+* •Wo(v 
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Thus, 


We have 


G n (x)= -\[Y n {x)-iJ n {x)] 






da; 

’ sin kxdx 
■y/x 2 — 1 

(XIV) 


TT v //\ . •! T/n f 00 COS kxdx , . 

Thus, 

f ® cos /tzdx 7r Tr /; N 

J, V3^T = ”^ oW 

(xv) 

Further, 


Thus, 



Vz 2 — 1 


(cos Arx 4- i sin Arz) da; 

= -^l (*) 

Fi® 


■yjx 2 — 1 
0 a; cos A;a;da; 

~~ 2 

® a; sin Arzdz 


(XVI) 

(XVII) 


TOTAL AERODYNAMIC FORCES AND MOMENTS 


TOTAL FORCE 

From equations (I) and (VIII) we obtain 
P = — p6 2 (V7ra + Trh — irbda — vl\ft — Tib ft) 

— 2-wpvbC | va 4- h 4- b^ 2 — 4- —Tuft ft 

+ b^T u , 3 1 (XVIII) 


TOTAL MOMENTS 

From equations (II) and (IX) we obtain similarly 

- P 6^{ -2T,-T 1 +T i (a-^vbdc+2T l3 b*& 

+ Wr, - r 4 r I0 ) - ±vbpT t T n - 3 

7T Z7T 7T 

— 2i — pvb 2r T\< 2 p | va 4- h 4- b (^2 — (i^d 

+\T w vfi+b±Tnb\ (XIX) 


From equations (III) and (X) 

M. = - pblir(~ - a)vba + 7 rfe 2 (g + a 2 )a 
+ (T 4 4- Ti 0 )v 2 ft 

+ (^l\-T 6 -(c-a)T 4 + ±T u ybp 

— ^ T 7 4- (c — a) b 2 ft — a7r6/tj 

4- 2pvb 2 ir^a 4- va + h- i- 

+ lT l0 v(i + b±T n ^ (XX) 

DIFFERENTIAL EQUATIONS OF MOTION 

Expressing the equilibrium of the moments about a 
of the entire airfoil, of the moments on the aileron 
about Cy and of the vertical forces, we obtain, respec- 
tively, the following three equations; 
a \ — I a d — Ipfi — b(c — a)Spfi — Seth ~ aCa 4- M a = 0 

ft 1 — I aft — I fid — b{c — d)aSfi — hSfi — (3Cfi 4- Mfi — 0 

h : - hM ~ dSet - 1 BSfi - hC h 4- P = 0 

Rearranged : 

a \ ala 4- ft ( Ifi + b (c — d)Sfi) 4- hS a + aC a — Al a = 0 

ft'. a(Ifi~h b(c~ d) Sfi) + ft Ifi 4 hSfi 4* ftCfi — Mfi = 0 

h : dS a 4- ftSfi 4- h,M+ KC h - P = 0 

The constants are defined as follows: 

p, mass of air per unit of volume. 

by half chord of wing. 

My mass of wing per unit of length. 

SaySfi, static moments of wing (in slugs-feet) per 
unit length of wing-aileron and aileron, 
respectively. The former is referred to 
the axis d ; the latter, to the hinge c. 
I a ,h, moments of inertia per unit length of 

wdng-aileron and aileron about a and c, 
respectively. 

C a y torsional stiffness of wing around d, cor- 

responding to unit length. 

Cf, torsional stiffness of aileron around c, cor- 

responding to unit length. 

C h) stiffness of w T ing in deflection, correspond- 

ing to unit length. 

DEFINITION OF PARAMETERS USED IN EQUATIONS 


Pb 2 

M 


the ratio of the mass of a cylinder of air of 
a diameter equal to the chord of the 
wing to the mass of the wing, both taken 
for equal length along span. 


408269 O 41 — 2 
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ra = ylm f ^ 1C rac ^ us gyration divided by b. 

S a 

x a = Alb’ center g ray ity distance of the wing 

from a. divided by b. 

ICa 

w a = 'W'j £L; the frequency of torsional vibration 
around a. 

r * = ^M&’ re( ^ uce< ^ radius of gyration of aileron 

divided by b, that is, the radius at 
which the entire mass of the airfoil 
would have to be concentrated to give 
the moment of inertia of the aileron Ip. 




Mb 

- lOl 

v 

/Cn 


reduced center of gravity distance from c. 

frequency of torsional vibration of aileron 
around c. 

frequency of wing in deflection. 


FINAL EQUATIONS IN NONDIMENSION AL FORM 


On introducing the quantities M a , Mp, and P, 
replacing T 9 and T i3 from page 5, and reducing to 
nondimensional form, we obtain the following system 
of equations: 


(A) 


(B) 


(C) 


«[^ 2 + «Q + a 2 )]+<i|-(Q-a)+a ] §- 2 +^[ 


P JL 

r l -f (c - a)xp k - (c - a) — 


1 x' t ] + ^l[“ 2 ^“G- a ) 7 *] 


+ @K L i 


Ji ( r 4 + ^+S(^- a «)J- a (a + |)^[^ + { + (i-a) A+ ^|/» + |^]-0 


«[V+ (c-a)^] +* (p- T. 


+ /3 


T \2 vC(k)Vv<x h /I 

K r — I — - 


a(*.-rt) + 45«+3(a»-Jr,«)-^2>i+X(l+«)J+A§J 


(!-a)d + ^|0+!^]-=O 


SOLUTION OF EQUATIONS 

As mentioned in the introduction, we shall only have 
to specify the conditions under which an unstable 
equilibrium may exist, no general solution being 
needed. We shall therefore introduce the variables at 
once as sine functions of the distance s or, in complex 

form with ^ as an auxiliary parameter, giving the 
ratio of the wave length to 2tt times the half chord b: 

ik T 

a = a 0 e ' 0 

and A = /t 0 e' 

where s is the distance from the airfoil to the first 
d s 

vortex element, = 0 , and and are phase angles 

of P and h with respect to a. 


Having introduced these quantities in our system of 
equations, we shall divide through by K - 

We observe that the velocity v is then contained in 
only one term of each equation. We shall consider 
this term containing v as the unknown parameter 9.X. 
To distinguish terms containing X we shall employ a 
bar; terms without bars do not contain X. 

We shall resort to the following notation, taking care 
to retain a perfectly cyclic arrangement. Let the 
letter A refer to the coefficients in the first equation 
not containing C(k ) or X, B to similar coefficients 
of the second equation, and C to those in the third 
equation. Let the first subscript a refer to the first 
variable a , the subscript j 8 to the second, and h to the 
third. Let the second subscripts 1, 2, 3 refer to the 
second derivative, the first derivative, and the argu- 
ment of each variable, respectively. A aX thus refers 
to the coefficient in the first equation associated with 
the second derivative of a and not containing C(Jc) or 
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X; C h 3 to the constant in the third equation attached to 
k, etc. These coefficients 4 are as follows: 


^= r ; 2 +G +a2 ) 

A a2 = (j-a) 

A a3 — 0 

p K 7 r \ K 7T / 

A »=l[- 2 P-(k~ a ) Tt ] 


Afa — (1 4 + 1 io) 



B « , = 0 


5 ) 


(-An) 


The solution of the instability problem as contained 
in the system of three equations A, B, and C is given 
by the vanishing of a third-order determinant of com- 
plex numbers representing the coefficients. The solu- 
tion of particular subcases of two degrees of freedom 
is given by the minors involving the particular co- 
efficients. We shall denote the case torsion-aileron 
(a, 0) as case 3, aileron-deflection (0, h) as case 2, and 
deflection-torsion ( h , a) as case 1. The determinant 
form of the solution is given in the major case and in 
the three possible subcases, respectively, by: 


D= 


- . V . 

Raa~\~H act) RaP~\~H uP) Rah~\~^lah 

Rba-\-Hba) Rb&-\riIbP t + 

R ca~\~H ca) Rcfi-\-Ucp) Rch~\~Hch 


-0 


and 


■jrj- I Raa~\~Haa) RaP~\~ H uP I ~ 

I Rba-\~Hba, Rb&-\-lIba I 
-rj _ | Rb&~\-H bp) Rbh J rHb)i I _ ~ 
\RcP~\~HcP) Rch~\~Hch\ 

-jTj I Rch~\~H Ch) Rca~\~H ca . „ 

M " l Ran+iU Ra. + ilJ 


Case 3 


Case 2 


Case 1 


L> _ r/3 L rp 


&P2 2^2 r R^ r R\\ 

B^^{T b -T t T w ) 


REAL EQUATIONS IMAGINARY EQUATIONS 


RaaRaf) 

RbaRbP 

- 

1 J-aaJ-aP 
1 IbalbP 

= 0 

1 RaaR ap 
1 1 bal bP 

+ ! halap 
1 1 RbaRbp 

= 0 Case 3 

RbpRbh 1 
RcpRcn 1 

~ 

l b&l bh j 
1 I cpl ch \ 

H 

RbpRbh 
I cpl ch 

1 , j Ibplbh j 

1 1 R cpRch 1 

= 0 Case 2 

RchRca 

RahRaa 

~ 

1 I chi ca 
\ 1 ahl-aa 

Uo 

j RchRca 

1 lahl aa 1 

. I chi ca I 

RahRaa 1 

=0 Case 1 


B hl — 0 

B/a = 0 

C.t-^-a (=A m ) 

C a 2 = 1 

Crt == o 

c °'= x i~l T ' < =B *> 

= —\t, 

Ch - 0 

C„ = ^ + i 

Ch2 — 0 
ChZ — o 

* The factor 1 or d i S no t included in these constants. See the expressions for 
the R’s and /’ s on next page. 


Note.— T erms with bars contain X; terms without bars do not contain X. 

The 9 quantities- R aa , R a? , etc., refer to the real parts 
and the 9 quantities I aaf I afit etc., to the imaginary 
parts of the coefficients of the 3 variables a, 0, and h 
in the 3 equations A, B , C on page 10. Denoting the 
coefficients of a, a , and a in the first equation by p, 
q , and r, 

which, separated in real and imaginary parts, gives 
the quantities R aa and I aa . Similarly, the remaining 
quantities R and I are obtained. They are all func- 
tions of fc or C(k). The terms with bars R aa , RbP f 
and R ch are seen to be the only ones containing the 
unknown X. The quantities il and X will be defined 
shortly. The quantities R and I are given in the 
following list: 
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A.= - A aI + X+ 1 + a) [(! - a)<?- 1 /] (1) 

' Ra)=— ^0l+p.A/S3+ + 2 jTioFJ (2) 

n ah =-A H + j2[a + ^G (3) 

( 4 ) 

, A S =-B ft +pB M H-O,X-Igi[r n e-2r 10 ^] (5) 
— Uni— (6) 


«^=-c sl -ii[r u <7-2r 10 |/] 

R C h= — Chi-\-SI) I X.—^2G 


(V) 

(8) 

(9) 


^ [ 2 ( a + i) I G ■ 0 ) F+ iM - a ° 2 ~1 
l[K 0+ i)( T " F+2 i ri » (? )- A *] 
H a+ i) F 


j _ A 

Iba k 

m 

_ 7r l 

( l 2~ a ) F+ Tc G \+ B °\ 

T _ 1 

hfi ~k 

7j2 / 
2tt 2 \ 

V 11 F+2|r 10 6f)+B w 

T — 4 "^ l2 Tf 




(id 

( 12 ) 

(13) 

(14) 

(15) 

(16) 


ics> = l l( TnF+2 l t '°°) + °»\ 

U = \2F 


(17) 

(18) 
(19) 


The solution as given by the three-row determinant 
shall he written explicitly in X. We are immediately 
able to put down for the general case a cubic equation 
in X with complex coefficients and can easily segregate 
the three subcases. The quantity D is as before the 
value of the determinant, but with the term containing 
X missing. The quantities M aa , M b0) and M ch are 
the minors of the elements in the diagonal squares 
a<x, b(3, and ch, respectively. They are expressed ex- 
plicitly in terms of R and I under the subcases treated 
in the following paragraphs. 


A au ~]r^l a X A /l0 


D= 


A bl 

A c< 


A b0 -\-Q 0 X 

A-cp 


Aah 

A bb 

A C h-\-^hX. 


=0 


where A aa =R aa +iI aa etc. 


Complex cubic equation in X: 


ilafytikX 3 + (Q a tt 0 A cn + Q 0 ft b A a(X + A b0 ) X 2 

+ (VaM aa +n fi M b ( i +to h M ch )X+D=Q (XXI) 

Case 3, (a, &) : 


W 2 + (n a A b0 +%A aa )X+M ch =Q 
Case 2, (0, h) : 

il 0 %X 2 -\- (® 0 A ch -\rtt b A b0 )X-{-M aa =O 
Case 1, (h, a) : 


VAX 2 + (VnAao+^A^X+M^O 


C a = / c o a r a V 1 /6r r to r y 

\o ) T r T ) K\vk) 


QaX~ 


QfiX- 


k 2 Mv 2 K 

C t 

k 2 Mv 2 K 

Cub 2 


_/w 0 r 0 \ 2 l /br r u T \ 2 
\w r r r / K\vkJ 

= ( V 1 ( br T o> T \ 2 

k 2 Mv 2 K \co T r T J k\ vk J 


and finally 


y 

k\ vk / 


(XXII) 

(XXIII) 

(XXIV) 


We are at liberty to introduce the reference param- 
eters w r and r T} and the convention adopted is: o> r is 
the last w in cyclic order in each of the subcases 3, 2, 
and 1. 


Then = ^ and i2 n+1 = l, thus for 

\co„ +1 r n+1 / ’ 

Case 3, and tt 0 = 1 

Case 2, and ft A =l 

Case 1, 8*=^"^ and il a = 1 


To treat the general case of three degrees of freedom 
(equation (XXI)), it is observed that the real part 
of the equation is of third degree while the imaginary 
part furnishes an equation of second degree. The 
problem is to find values of X satisfying both equa- 
tions. We shall adopt the following procedure: Plot 

graphically X against | for both equations. The points 
of intersection are the solutions. We are only con- 
cerned with positive values of ^ and positive values of 

X. Observe that we do not have to solve for k , but 
may reverse the process by choosing a number of 
values of k and solve for X. The plotting of X 

against ^ for the second-degree equation is simple 

enough, whereas the task of course is somewhat more 
laborious for the third-degree equation. However, 
the general case is of less practical importance than 
are the three subcases. The equation simplifies con- 
siderably, becoming of second degree in X. 
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We shall now proceed to consider these three sub- 
cases. By virtue of the cyclic arrangement, we need 
only consider the first case ( a , 0). The complex 
quadratic equations (XXII) ~ (XXIV) all resolve 
themselves into two independent statements, which 
we shall for convenience denote “Imaginary equa- 
tion” and “Real equation”, the former being of first 
and the latter of second degree in X. All constants 
are to be resolved into their real and imaginary parts, 
denoted by an upper index R or /, respectively. 

Let M aa =M R a<x -\-iM I aa and let similar expressions 
denote M b p and M ch 

Case 3, (a,p). Separating equation (XXII) we obtain. 
(1) Imaginary equation: 

(i2 «/ b p -J- ttplaa) X T M 1 ch = 0 

y 

ti a Ibp-\-topIa<x 

(21 Real equation: 

ftaftflA 2 -)- (&ctRbp-\-ttpRaa)X-\-M R c }, = 0 

Eliminating X we get 

tiA(M T ch ) 2 — (Qjl b p+QpR aa ) (SlJbP+QpIaaW 1 * 
+M R ch (sij bfi +sipi aa y = o 

By the convention adopted we have in this case: 

°- — Cs)’CS)** and Q * =1 

Arranging the equation in powers of ft a we have: 

— A t 1 c h(R aal 6/9 + 1 actRbp) ~\~2M R ch l aa I 6^] 

+ [-M I ch R aa I aa +M R ch I aa 2 ] = 0 
But we have 


(M 1 c*) 2 — M 1 C h (RaalbP 4“ IaaRbfi) 

— Af^ ch\RaaI bp RapI ba~\~ R b pl aa R balaP Raa I bp Rbp I a«] 
= — A/ 7 ch(RapI 6«+ 1 apRba) 

Finally, the equation for Case 3 (a, /3) becomes: 

0 a 2 (M R C J bp 2 — A i 1 chR bplbp) + [ — A t 1 ch (R a pl ba + IapR ba ) 

+2M R c JaJbp]+M R c J a J-M' ch R a J aa =0 (XXV) 

where 

Af ch == RaaRbP RapRba I aal bp~\~ Iaplba 

Af ch z= Raaf bP Rapt 6a"f laaRbP lapRba 

The remaining cases may be obtained by cyclic 
rearrangement: 


Case 2, (jS,^) co r =co A ^ = (cir) 

Bj(A l R a0 Jch MiJtcJen) -j- ft#! MURbJcp + lbbRcp) 
+2M R a JbpI ch ]±M R J b ?-M{MbpIbP=Q (XXVI) 
where M R aa =R b pRr,h— R bh R c p—I b pIch+IbhI c p 

^1 aa ==: R bpl ch Rbh^cP I bpRch I bhRcP 

Case 1, (A, «) « r =a>„ a h =(^j A fi„=i 

W h (M* h pIaS A I r b pR aa I aot ) -J- ft*[ — Mlp (R ca I ah + 1 caRah) 
+2M%I c J aa ]+M%IJ-mpRchI C h=0(XXVII) 


where A I b p — RchRaa RcaRah I chi aa~\~ I cal ah 

M b p=R C hIaa R cal ah~\~ I cb R aa I caRah 

Equations (XXV), (XXVI), and (XXVII) thus 
give the solutions of the cases: torsion-aileron , aileron- 
deflection, and deflection-torsion , respectively. The 
quantity ft may immediately be plotted against 

^ for any value of the independent parameters. 

The coefficients in the equations are all given in terms 
of R and /, which quantities have been defined above. 
Routine calculations and graphs giving ft against 

^ are contained in Appendix I and Appendix II. 

Knowing related values of ft and p X is immediately 

expressed as a function of ft by means of the first- 
degree equation. The definition of X and ft for each 
subcase is given above. The cyclic arrangement of 
all quantities is very convenient as it permits identical 
treatment of the three subcases. 

It shall finally be repeated that the above solutions 
represent the border case of unstable equilibrium. 
The plot of X against ft gives a boundary curve between 
the stable and the unstable regions in the Aft plane. 

It is preferable, however, to plot the quantity p 


instead of X, since this quantity is proportional to the 
square of the flutter speed. The stable area can easily 
be identified by inspection as it will contain the axis 


J. 1 
k 2 


-=£=0, if the combination is stable for zero velocity. 


Langley Memorial Aeronautical Laboratory, 
National Advisory Committee for Aeronautics, 
Langley Field, Va., May 2, 1934 • 
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APPENDIX I 


PROCEDURE IN SOLVING NUMERICAL EXAMPLES 


(1) Determine the R’s and I’s } nine of each for a 
major case of three degrees of freedom, or those per- 
taining to a particular subcase, 4 R’s and 4 Fs. Refer 
to the following for the R’s and Fs involved in each 
case: 

The numerals 1 to 9 and 11 to 19 are used for con- 
venience. 

(Major case) Three 
degrees of freedom 


1 

Raa 

laa 

11 

2 

R a p 

IaP 

12 

3 

Rah 

I ah 

13 

4 

Rba 

I ba 

14 

5 

Rbp 

IbP 

15 

6 

Rbh 

Ibh 

16 

7 

Rea 

lea 

17 

8 

Rep 

IcP 

18 

9 

Rch 

Ich 

19 


(Case 3) Torsional- 
aileron ( a , /?) 

1 

R aa 

I aa 

11 

2 

Rap 

lap 

12 

4 

Rba 

I ba 

14 

5 

RbP 

IbP 

15 


(Case 2) Aileron- 
deflection (/?, h) 

5 

Rbp 

IbP 

15 

6 

Rbh 

Ibh 

16 

8 

Rcp 

IcP 

18 

9 

Rch 

Ich 

19 


(Case 1) Deflection- 
torsion Qi, a) 

7 

Rea 

lea 

17 

9 

Rch 

Ich 

19 

1 

Raa 

I aa 

11 

3 

Rah 

I ah 

13 


It has been found convenient to split the R’s in two 
parts R=R'+R", the former being independent of 

the argument ^ . The quantities I and R" are func- 


tions of the two independent parameters a and c only. 6 
The formulas are given in the following list. 


BV=j2(a + |) ((!-«)<?-£) 

(1) 

R"^\l\(T l+ r ,„) !+(«+ 1 t k f) 

| (2) 

R " ah = lc 2 { a ' + ^) G 

(3) 


(4) 

w ! Y (?■„(?- f (7’ 5 - t, r, 0 ) 

| (5) 

nrr 1 ^12 p 

(6) 

kil"* 

1 

1 

r-qc^ 

CM 

1 

e 

03 

(7) 

to 

•» 

II 

1 

»-* 

^ 1 - 

* 

1 

to 

(8) 

R" ch— 2 G 

(9) 

/.„=-2(a + i)((i-a)F+itf) + |-a 

(11) 

/ oS =-^(a+|)(2’nF + |r 10 G)+ 2 P 

(12) 

+ 6 _ “) r< ) 


7 0 „=-2(a + |)F 

(13) 

Where p = -~(l-c 2 ) 3/2 

(14) 


h e = ^\T»(T n F+\T w G)-T t T n 

I - R —F 
Ith ~ v 1 


(15) 

(16) 


ft 

II 

to 

j- a ) F +H +1 


'T n F+\T m G)-r 

Ln = 2F 



(17) 

(18) 
(19) 


s The quantities I given in the appendix and used in the following calculations 
are seen to differ from the I’s given in the body of the paper by the factor ~ k • It 
may be noticed that this factor drops out in the first-degree equations. 


14 
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Choosing certain values of a and c and employing 
the values of the T’s given by the formulas of the report 
(p. 5) or in table I and also using the values of F and 
G (formulas (XII) and (XIII)) or table II, we evaluate 

the quantities I and K" for a certain number of ^ 

values. The results of this evaluation are given in 
tables III and IV, which have been worked out for 
a = 0,— 0.2, and — 0.4, and for c=0.5 and c = 0. The 

range of ^is from 0 to 40. These tables, save the work 

of calculating the I’s and R"’s for almost all cases of 
practical importance. Interpolation may be used for 
intermediate values. This leaves the quantities R' to 

be determined. These, being independent of p are as 

a result easy to obtain. Their values, using the same 
system of numbers for identification, and referring to 
the definition of the original independent variables on 
pages 9 and 10, are given as follows: 

^--T-(l +a2 ) <» 

«'-i— +5 + («-*) T (2) 

7 + a (3) 


Case 3 


Find products 1.5 2.4 11.15 12.14 

Thenil4 B c# =1.5— 2.4— p(l 1.15— 12.14) 

Find products 1.15 2.14 11.5 12.4 

Then M' c »=1.15— 2.14+11.5— 12.4 
and a=M B e# (15) 2 — M 7 e ,(5.15) 

6 = — M' c# (2. 14 + 12.4) + 2M fi e *(l 1 . 1 5) 
c=M R C h(liy — M r a (l.l 1) Find a. 

Solution :> 

A M r ch 

Similarly 


Case 2 


and 


5.9 6.8 15.19 16.18 

M R „=5.9— 6.8 — p(15.19— 16.18) 

5.19 6.18 15.9 16.8 

M , „„=5.19-6.18+15.9-16.8 
a=Af* 0 „(19) 2 — M 7 „„(9.19) 
6=-2l/'„„(6.18+16.8) 

+ 2M B „„(6.18+16.8) 
c=M R aa {l5) 2 —M I aa (5.15) Findft^ 

1 ft*(19) + 15 

A M\ a 

Case 1 


i?' &a ==same as R' a p (4) 

(5) 

--* + +, (6) 

k 7 r 

R ' cu = same as R' ah (7) 

R' C fi = same as R' hh (8) 

1 (9) 

Because of the symmetrical arrangement in the 

determinant, the 9 quantities are seen to reduce to 

6 quantities to be calculated. It is very fortunate, 

indeed, that all the remaining variables segregate them- 
selves in the 6 values of R' which are independent of p 

while the more complicated I and R" are functions 
solely of c and a. In order to solve any problem it is 
therefore only necessary to refer to tables III and IV 
and then to calculate the 6 values of R'. 

The quantities (1) to (9) and (11) to (19) thus 

having been determined, the plot of ft against p which 

constitutes our method of solution, is obtained by 
solving the equation aft 2 +5ft+c=0. The constants 
a, 6, and c are obtained automatically by computation 
according to the following scheme: 


9.1 

M R bP = 9. 1-7.3 
9.11 


7.3 19.11 17.13 

~p(19. 11-17. 13) 

7.13 19.1 17.3 


MV=9.11 — 7.13 + 19.1 — 17.3 
a=M R b0 (ny-M I bP (l.ll) 
b= — M / 6/S (7.13 + 17.3) + 2M^(19.11) 
c = M R bP ( 1 9) 2 — M 1 b p (9.19)* Find ft„ 


1 _P-*(11) + 19 
X ~ M 2 bfi 


ft a is defined as for case 3 ; 

\uprpj 

Up is defined as for case 2 ; and 

ft* is defined as ( ) f or case l . 

\co«r«/ 

The quantity is by definition. 

Since both ft and are calculated for each value of 

jp w T e may plot p directly as a function of ft. This 

quantity, which is proportional to the square of the 
flutter speed, represents the solution. 

W e shall sometimes use the square root of the above 

quantity, viz, -r -\l -L = * and will denote this 

K \ A b<jo T r r 


42 


16 


REPORT NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS 


quantity by F, which we shall term the “flutter factor.” 
The flutter velocity is consequently obtained as 


v=F b -^ 

Vk 

Since F is nondimensional, the quantity must 

V* 

obviously be a velocity. It is useful to establish the 
significance of this velocity, with reference to which 

the flutter speed, so to speak, is measured. Observing 
12 

that k = and that the stiffness in case 1 is given by 


(j) a 


I C ° 

V Mbh* 


this reference velocity may be written: 


■yJK 0 \ 7T p 


or 


T P V R *b 2 =C a 


The velocity v R is thus the velocity at which the total 

force on the airfoil irpv R 2 2b attacking with an arm ^ 

equals the torsional stiffness C a of the wing. This 
statement means, in case 1, that the reference velocity 
used is equal to the “divergence” velocity obtained 
with the torsional axis in the middle of the chord. This 
velocity is considerably smaller than the usual diver- 
gence velocity, which may be expressed as 


1 

Vd — Vr J 

2 + o 

where a ranges from 0 to — We may thus express 
the flutter velocity as 

v F — v r F 


In case 3 the reference velocity has a similar signifi- 
cance, that is, it is the velocity at which the entire lift of 

the airfoil attacking with a leverage ~ b equals numeri- 
cally the torsional stiffness Cp of the aileron or movable 
tail surface. 

In case 2, no suitable or useful significance of the 
reference velocity is available. 

TABLE I.— VALUES OF T 



c=i ! 

c=Yz 

c=0 

H 

c— — 1 

T,_._ 

0 

-0. 1259 

-0. 6667 

-1. 6967 

-3. 1416 

Tt 

0 

-0. 2103 

-1.5707 

-4. 8356 

-9. 8697 

Ti 

0 

— . 05313 

-.8084 

-3. 8375 

-11. 1034 

Ti 

0 

-.6142 

-1.5708 

-1.6614 

-a. i4i6 

Ti — 

0 

9398 

-3. 4674 

-6. 9503 

-9. 8697 

T t . 

0 

-0. 2103 

-1. 5707 

-4. 8356 

-9. 8697 

Ti 

0 

.0132 

-. 1964 

-1. 1913 

-3. 5343 

Ti 

0 

.0903 

-. 3333 

-1.4805 

-3. 1416 

Tu 

0 

1. 9132 

2. 5708 

2. 9604 

3. 1416 

T n - 

0 

1. 2990 

3. 5708 

6. 3538 

9. 4248 

Tit 

0 

. 07066 

. 4292 

1. 2990 

3. 1416 


TABLE II.— TABLE OF THE BESSEL FUNCTIONS J 0 , 
To, Yi AND THE FUNCTIONS F AND G 


m 


_./,(■/, + Yq) + Y,( Y 1 -Jo ) 
{ji+Yt)*+(Yi-Jt)* ' 


-<?(*)- 


Wi+FoWitri-Jo) 
(Ji + Vo)2+(V|-Jo)* 


: k 

1 

k 

Jo 

Ji 

Vo 

Yi 

F 

-G 

00 

0 





0. 5000 

0 

10 

Ho 

-0. 2459 

0. 0435 

0. 0557 

0.2490 

..5006 

0.0126 

6 

H 

. 1506 

-. 2767 

-.2882 

-. 1750 

.5018 

.0207 

4 

H 

-. 3972 

-. 0660 

-. 0170 

. 3979 

.5037 

. 0305 

2 

tt 

. 2239 

.5767 

.5104 

-. 1071 

.5129 

. 0577 

1 

l 

. 7652 

.4401 

.0882 

-. 7813 

.5395 

.1003 

.8 

1H 

.8463 

.3688 

-.0868 

-.9780 

. 5541 

.1165 

.6 

1H 

. 9120 

.2867 

-.3085 

-1.2604 

.5788 

. 1378 

.5 

2 

. 9385 

.2423 

-. 4444 

-1.4714 

.6030 

.151 

.4 

2 tt 

. 9604 

. 1960 

-. 6060 

-1. 7808 

. 6245 

. 166 

.3 

3tt 

.9776 

.1483 

-. 8072 

—2. 2929 

.6650 

. 180 

.2 

5 

. 99(H) 

. 0995 

-1.0810 

-3. 3235 

. 7276 

. 1886 

.1 

10 

. 9975 

.0499 

-1.5342 

-7.0317 

.8457 

. 1626 

.05 

20 





. 911 

. 132 

. 025 

40 





. 965 

.090 

0 






1.000 

0 
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TABLE III.— VALUES OF R 



1 

k 

c 

a 

0 

Ho 

H 

H 

H 

■ 

1H 

m 

2 

2H , 

3H 

5 

10 

20 

40 



0 

0 

-0.00564 

-0.01566 

-0. 03529 

-0. 14265 

-0. 58965 

-0.93656 

-1. 72330 

-2. 56300 

-4.11000 

-7.68720 

-18. 66150 

-85.38300 

-365.72000 

-1,528. 2000 

R"aa 

0) 

-.2 

0 

-.00353 

-.00981 

-.02208 

-.08905 

-. 36586 

-.58061 

-1.08158 

-1.57400 

-2. 51580 

-4.68430 

-11.31010 

-51. 42490 

-219.74900 

-917. 3520 



-.4 

0 

-. 00123 

-. 00341 

-. 00767 

-.03084 

-. 12595 

-. 19936 

-. 36305 

-. 53676 

-. 85520 

-1. 58540 

-3. 80774 

-17. 20670 

-73.35520 

-305. 9280 



0 

0 

00163 

-.00452 

-.01020 

-. 04175 

-. 18016 

-. 29384 

-. 56223 

-.87212 

-1. 43983 

-2. 84988 

-7. 46300 

-38. 29650 

-172.36360 

-741. 7972 


0 

-.2 

0 

.00030 

.00083 

.00184 

.00679 

. 01922 

. 02266 

.01629 

-.01400 

-.06803 

-. 29517 

-1.29480 

-10. 24590 

-52. 49020 

-241. 3664 

P" * 


-.4 

0 

.00222 

.00617 

.01388 

.05531 

. 21861 

. 33914 

. 59499 

. 84414 

1. 30365 

2. 25914 

4.87340 

17. 80470 

67. 38320 

259. 0648 

«p 


0 

0 

.00083 

.00229 

.00510 

. 01932 

.06419 

.08876 

. 12176 

. 12260 

.12205 

-.02900 

-.93535 

-10. 48970 

-59. 16180 

-268. 7236 


0.5 

-.2 

0 

.00214 

.00595 

.01336 

.05278 

.20325 

.31065 

. 53062 

.73222 

1. 10233 

1.81135 

3. 55230 

10. 14740 

31. 49620 

101. 6340 



-.4 

0 

.00347 

.00965 

.02170 

.08656 

. 34361 

.53463 

.94336 

1. 34762 

2. 09190 

3. 66913 

8.08235 

30. 97980 

120.89760 

475.2592 



0 

0 

-.00125 

-.00345 

-. 00763 

-.02890 

-. 10030 

-. 14560 

-. 22470 

-.30200 

-.41500 

-.60000 

-.94300 

-1.62600 

-2. 64000 

-3.6000 

R"ab 

(0 

-.2 

0 

-. 00075 

-. 00207 

-.00426 

-. 01734 

-. 06018 

-. 08736 

-. 13482 

-. 18120 

-. 24900 

-.36000 

-. 56580 

-.97560 

-1.58400 

-2. 1600 



-.4 

0 

-.00201 

-. 00334 

-. 00502 

-. 01003 

-.02006 

-.02508 

-. 03236 

-. 04012 

-.05015 

-. 06683 

-. 10030 

-. 20060 

-.40120 

-. 8024 



0 

0 

.00077 

.00214 

.00482 

. 01949 

.08055 

.12821 

.23541 

.35010 

. 56143 

1.05008 

2.54920 

11.66330 

49.95700 

208. 7520 


0 

-.2 

0 

.00080 

.00223 

.00503 

.02027 

.08329 

. 13219 

. 24169 

. 35836 

. 57276 

1. 06650 

2. 57490 

11.70770 

50.03000 

208. 8500 

R"ba 


-.4 

0 

.00084 

.00233 

.00523 

.02106 

.08603 

. 13616 

. 24796 

. 36661 

.58410 

1. 08286 

2. 60069 

11. 75220 

50. 10160 

208. 9490 


0 

0 

.00013 

.00035 

.00079 

.00321 

.01327 

.02112 

.03878 

.05767 

.09248 

.17296 

. 4 1988 

1.92110 

8. 22870 

34. 3850 


0.5 

-.2 

0 

.00013 

.00037 

.00083 

.00334 

. 01372 

.02177 

. 03981 

, .05903 

.09434 

. 17566 

. 42413 

1. 92840 

8. 24060 

34. 4007 



-.4 

0 

.00014 

.00038 

.00086 

.00347 

.01417 

. 02243 

.04084 

.06039 

. 09621 

. 17836 

.42837 

1.93575 

8. 25246 

34. 4169 

R " bp 

0 


0 

.00124 

.00343 

.00772 

.03101 

.12642 

. 19830 

. 35807 

.52400 

.82930 

1.51630 

3. 54970 

15. 35120 

64. 02240 

263. 2340 


.5 

( s ) 

0 

.00031 

.00087 

.00196 

.00785 

.03170 

. 04980 

.08935 

.13000 

.20440 

. 36940 

. 84970 

3. 55050 

14. 56740 

59. 3188 

D//. . 

0 


0 

.00017 

.00047 

.00104 

.00394 

. 01370 

.01989 

. 03177 

.04125 

.05669 

. 08196 

.12881 

. 22211! . 36062 

.4918 

n bh 

.5 

( 3 ) 

0 

.00003 

.00008 

.00016 

.00065 

. 00226 

.00328 

.00506 

.00680 

.00934 

.01350 

.02122 

. 03659 

.05940 

.0810 



0 

0 

.01128 

.03132 

.07058 

.28530 

1. 17930 

1. 87710 

3. 44670 

5. 12600 

8. 22000 

15. 37450 

| 37.32300 

170. 76600 

731. 44000 

3, 056. 4000 

R"ca 

<>) 

-.2 

0 

. 01178 

.03270 

. 07362 

. 29684 

1. 21954 

1. 93540 

3. 53860 

5. 24680 

8. 38600 

15. 61440 

37.70020 

171.41640 

732. 49600 

3, 057. 8400 



-.4 

0 

.01228 

. 03408 

.07668 

.30838 

1. 25950 

1. 99360 

3. 63050 

5. 36760 

8. 55200 

15. 85440 

38. 07740 

172. 06680 

733. 55200 

3, 059. 2800 

R " A 

0 

(*) 

0 

-.00963 

-. 02673 

-. 06018 

-. 24266 

-1.00561 

-1.58246 

-2. 89371 

-4.29100 

-6. 85898 

-15.49965 

-30. 84330 

-140.263701 -599.41300 

-2, 502. 3470 

n ep 

.5 


0 

.00660 

. 01840 

. 04150 

. 16810 

. 69850 

1. 11453 

2. 05320 

3. 06224 

4. 92530 

9. 24438 

22. 54400 

103. 67300 

444. 86400 

1, 881. 9900 

R"cb 

(') 

<*) 

0 

.00250 

.00690 

.01420 

. 05780 

.20060 

.29120 

. 44940| . 60400 

.83000 

1.20000 

! 1.88600 

3. 25200 

5. 28000 

7.2000 


1 Independent of c, * Independent of o. 


TABLE IV.— VALUES OF I 


1 

k 

0 

Ho 


Vi 

H 

J 

m 

m 

2 

2H 

3H 

5 

10 

20 

40 


C 

a 







0 

0. 25000 

0. 25096 

0. 25255 

0. 25578 

0.27240 

0. 33055 

0. 36855 

0.44030 

0.50050 

0. 60275 

0. 76750 

1.07920 

1. 70320 

2. 68450 

3. 61750 

/a. 

0) 

-.2 

.49000 

. 49050 

.49131 

. 49302 

.50189 

.53359 

. 55464 

. 69472 

. 62794 

.68671 

.78070 

.96021 

1. 32040 

1.90140 

2.45470 


-.4 

.81000 

.81014 

. 81037 

.81086 

.81145 

.82395 

.82938 

.84176 

. 85186 

. 87059 

.90030 

.95763 

1. 07300 

1. 26400 

1.44630 



0 

. 17805 

. 17874 

. 17985 

. 18219 

. 19433 

.23768 

.26645 

. 32132 

.36664 

.44690 

. 57526 

.82035 

1.31213 

2. 10476 

2. 85963 


0 

-.2 

.39170 

. 39212 

. 39278 

.39418 

. 40147 

. 42748 

. 44474 

.47761 

.50485 

.55300 

. 63002 

.77708 

1. 07215 

1. 54773 

2. 00065 

U 


-.4 

.60531 

.60545 

. 60567 

.60614 

. 60857 

. 61724 

.62299 

.63395 

.64303 

.65908 

.68475 

. 73377 

.82313 

. 99065 

1. 14163 


0 

. 13252 

. 13317 

. 13425 

.13640 

. 14742 

. 18544 

.20914 

.25611 

.29514 

. 35951 

. 46379 

. 65973 

1. 05124 

1. 65524 

2.22869 


0.5 

-.2 

.21297 

. 21336 

. 21401 

. 21530 

.22191 

. 24-172 

.25894 

.28712 

.31054 

. 34916 

.41173 

.52929 

.76420 

1. 12651 

1.47067 



-.4 

.29342 

.29354 

.29376 

.29419 

.29640 

.30-1.' 

.30891 

.31813 

.32594 

.33881 

. 35966 

.39884 

.47714 

. 59792 

.71260 



0 

-.50000 

-.50060 

-.50180 

-. 50370 

-.51290 

-. 53950 

-. 55410 

-. 57880 

-.60300 

-. 62450 

-. 66500 

-. 72760 

-. 84570 

-.94100 

-. 96500 

lab 

0) 

-.2 

-.30000 

-.30036 

-.30108 

- 30222 

-. 30774 

-.32370 

— . 33246 

-. 34728 

-. 36180 

-. 37470 

-. 39900 

-. 43656 

-. 50752 

-. 56460 

-. 57900 


-.4 

-.10000 

-.10012 

-. 10036 

-. 10074 

-. 10258 

-. 10790 

-. 11082 

-. 11576 

-.12060 

-.12490 

-. 13300 

-. 14552 

-. 16914 

-. 18220 

-. 19300 



0 

.39023 

.39010 

.38988 

. 38944 

.38717 

.37923 

.37404 

. 36424 

.35601 

.34204 

.31954 

.27696 

. 19172 

.05766 

-. 06980 


0 

-.2 

. 40389 

.40378 

.40359 

. 40320 

.40119 

. 39397 

. 38918 

.38005 

.37249 

.35911 

. 33771 

.29683 

.21469 

.08255 

-.04344 

r. 


-.4 

. 41755 

. 41746 

. 41730 

. 41697 

.41520 

.40871 

.40432 

.39586 

.38896 

.37617 

.35598 

.31671 

.23793 

.10744 

-. 01707 

i&a 


0 

.07438 

. 07435 

.07433 

.07424 

.07387 

.07256 

. 07171 

.07009 

.06874 

.06644 

.06273 

. 05572 

.04168 

.01960 

-.00327 


0.5 

-.2 

.07663 

. 07661 

.07658 

. 07651 

. 07618 

.07499 

.07420 

.07270 

. 07145 

.06925 

.06572 

.05899 

.04548 

.02370 

. 00295 



-.4 

. 07887 

. 07885 

. 07882 

. 07867 

. 07848 

. 07741 

.07668 

.07529 

. 07416 

.07205 

.06871 

.06226 

.04928 

.02779 

.00728 


0 

(*) 

.32297 

.32288 

. 32273 

.32241 

.32075 

.31483 

.31090 

. 30342 

. 29721 

.28625 

.26872 

.23524 

.16806 

. 05979 

-.04333 

.5 

.04270 

.04270 

. 04270 

.04270 

.04240 

.04150 

.04095 

.03930 

.03904 

.03760 

.03386 

.03080 

.02200 

.00845 

-.00470 

Ibb 

0 

(*) 

.06830 

.06840 

.06850 

.06880 

. 07010 

. 07370 

.07570 

.07910 

.08240 

. 08530 

.09080 

. 09940 

.11550 

.12440 

.13180 

.5 

. 01125 

. 01126 

.01129 

.01133 

.01154 

. 01214 

.01247 

.01302 

.01357 

. 01405 

. 01496 

.01637 

.01903 

.02117 

. 02171 



0 

1.50000 

1.49808 

1.49490 

1.48844 

1.45520 

1. 33890 

1.26290 

1. 11940 

.99900 

.74950 

.46500 

-.15840 

-1.40630 

-3. 36900 

-5. 23500 

lea 

(') 

-.2 

1.70000 

1.69832 

1. 69562 

1. 68992 

1. 66036 

1. 55470 

1. 48454 

1. 35092 

1. 24020 

1. 04430 

.73100 

.13264 

-1.06802 

-3. 00460 

-4. 84900 


-.4 

1.90000 

1. 89856 

1. 89634 

1. 89140 

1. 88552 

1. 77050 

1. 70618 

1. 58244 

1.48410 

1.31410 

.99700 

.42370 

-.72974 

-2. 64020 

-4. 46300 

I'fi 

0 

(’) 

1.06830 

1.06690 

1. 06470 

1.06000 

1.03580 

. 94900 

.89150 

.78190 

.69110 

.52840 

. 27380 

-. 21640 

-1.20010 

-2.78550 

-4. 29530 

.5 

. 40220 

.40100 

.39880 

.39450 

. 37240 

.29640 

.24640 

. 15510 

.07610 

-. 05180 

-.26030 

-.65220 

-1. 43520 

-2. 64380 

-3. 79010 

Ieb 

(>) 

(*) 

1.00000 

1.00120 

1.00360 

1. 00740 

1.02580 

1. 07900 

1. 10820 

1. 15760 

1. 20600 

1. 24900 

1.33000 

1. 45520 

1. 69140 

1.82200 

1. 93000 


1 Independent of e. J Independent of a. 
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APPENDIX II 

NUMERICAL CALCULATIONS 


A number of routine examples have been worked out | 
to illustrate typical results. A “ standard” case has 
been chosen, represented by the following constants: 

*=0.1, c=0.5, a=— 0.4, x a =0.2, 

rj= 0.25, Z/> = go’ V = i60 
w«, w/s , w* variable. 

We will show the results of a numerical computation 
of the three possible subcases in succession. 



Case 3, Torsion-aileron (a,/3): Figure 5 shows the 

against ^ relation and figure 6 the final curve 

F=k( — Yagainst S2«=(^Y=40^Y 
\o)prpbJ b \upri}/ \«/s/ 



a* 

Figure 6— Case 3, Torsion-aileron ( a , ff ) : Standard case. Showing flutter factor 
F against ll a . 

Case 2, Aileron-flexure (0, h ): Figure 7 shows the 
against relation 6 * and figure 8 the final curve 

against = 

6 It is realized that considerable care must be exercised to get these curves reason- 

ably accurate. 

18 


The heavy line shows the standard case, while the 
remaining curves, show the effect of a change in the 

value of xp to and — tq^ 

Case 1, Flexure-torsion ( h , a): Figure 9 shows again 



Figure 7.— Case 2, Aileron-deflection (fl, ft) : (a) Standard case, (b), (c), (d) indicate 
dependency on xp. Case (d), xp= —0.004, reduces to a point. 

the against ^ relation and figure 10 the final result 

* feb)’ "« ainst n * = fe) 8 “ 4 © 2 

Case 1, which is of importance in the propeller theory, 
has been treated in more detail. The quantity F shown 

in the figures is y* — yy 

Figure 11 shows the dependency on — = ^ ; 

figure 12 shows the dependency on the location of the 
axis a; figure 13 shows the dependency on the radius of 
gyration r a —r; and figure 14 shows the dependency 
on the location of the center of gravity x, for three 
different combinations of constants. 

EXPERIMENTAL RESULTS 

Detailed discussion of the experimental work will not 
be given in this paper, but shall be reserved for a later 
report. The experiments given in the following are 
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Wing B, wood, with flap, and the constants: 

c= 0.5, a=— 0.4, x«=0.192, r* 2 =0.178, 

£0=0.019, ?*0 2 =O.OO79, and co a kept constant 

= 17.6X2tt 

The results for wing A, case 1, are given in figure 15; 
and those for wing B, cases 2 and 3, are given in figures 
16 and 17, respectively. The abscissas are the fre- 
quency ratios and the ordinates are the velocities in 
cm/sec. Compared with the theoretical results calcu- 
lated for the three test cases, there is an almost perfect 
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Figure 12.— Case 1, Flexure-torsion (ft, a) : Showing dependency of F on location 
of axis of rotation a. Airfoil with r=4-lx =o.2;/t =h; — = a variable. 

2 4 ui 6 

agreement in case 1 (fig. 15). Not only is the minimum 
velocity found near the same frequency ratio, but the 
experimental and theoretical values are, furthermore, 
very nearly alike. Very important is also the fact that 
the peculiar shape of the response curve in case 2, pre- 
dicted by the theory, repeats itself experimentally. 
The theory predicts a range of instabilities extending 
from a small value of the velocity to a definite upper 
limit. It was very gratifying to observe that the upper 
branch of the curve not only existed but that it was 
remarkably definite. A small increase in speed near 
this upper limit would suffice to change the condition 
from violent flutter to complete rest, no range of transi- 
tion being observed. The experimental cases 2 and 3 
are compared with theoretical results given by the 
dotted fines in both figures (figs. 16 and 17). 


The conclusion from the experiments is briefly that 
the general shapes of the predicted response curves re- 
2.00 r 
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Figure 13.— Case 1, Flexure-torsion (A, a) : Showing dependency of F on the radius 
of gyration r«=r. 

A, airfoil with a= -0.4; «=-Ii z=0.2> r variable. 

B, airfoil with a= —0.4; *=4-’ i=0.25 —=1.00; r variable. 



A, airfoil with r 

B, airfoil with r 

C, airfoil with r= 


location of the center of gravity. 
1 . 

a=— u.% — =" ~77 

2 400 on 6 

1 • „ . 1 • 


1 . 


x variable. 

4-5 <z=— 0.4; “ I ’* variable. 

2 4 o>2 6 


1 . 


a = -0.4; <c=i 7 ^ ; — “1; x variable. 


peat themselves satisfactorily. Next, that the influ- 
ence of the internal friction 7 obviously is quite appreci- 


^ This matter is the subject of a paper now in preparation. 
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21 


able in case 3. This could have been expected since 
the predicted velocities and thus also the air forces on 
the aileron are very low, and no steps were taken to 
eliminate the friction in the hinge. The outline of the 
stable region is rather vague, and the wing is subject 



Figure 15.— Case 1. Wing A. Theoretical and experimental curves giving flutter 
velocity v against frequency ratio — • Deflection-torsion. 


to temporary vibrations at much lower speeds than 
that at which the violent flutter starts. The above 
experiments are seen to refer to cases of exaggerated 
unbalance, and therefore of low flutter speeds. It is 
evident that the internal friction is less important at 
larger velocities. The friction does in all cases increase 
the speed at which flutter starts. 



w 9 /u) h 


Figure 16.— Case 2. Wing B. Theoretical and experimental curves giving flutter 
velocity v against frequency ratio — • Aileron-deflection 08, A). 
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Figure 17.— Case 3. Theoretical curve giving flutter velocity against the fre- 
quency ratio — • The experimental unstable area is indefinite due to the im- 

<A>0 

portance of internal friction at very small velocities. Torsion-aileron (a, /9). 
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APPENDIX III 


EVALUATION OF 


r- 


(x— Xi) 2 +(y— y,) 2 
(x— x,) 2 +(y+y,) 2 


dxi 


*l) 2 +(l/ + 2/l) 2 Je 


2 y 


r 


Xidx t 
V\ (x—xt) 


, C l XjdXj (' dx t 

Jc ^l—x{ 2 (x — Xi) J VI— ^l 2 

J * ^ 3 . 

7 r- }z 2 [Putting X \ = cos 0] 

(Zi— x) VI— *i 

x 1 —x cos 0+ -y 1 — x 2 sin 0 j cnsg=l 

Vl— x 2 ° g cosO— x Los e= c 

i — cx-\- yi — x 2 yi — c 2 


= COS" 1 C + 


Vi-r 


;log- 


= COS —1 C “I - - 


,l°g: 


c—x 

c—x 


Vl — X 2 ® 1— cx— Vl— X 2 ^j\— c 2 
(p ~ = — %c log (1— cx— Vl— x 2 Vl — c 2 )+2c log (x— c) 
— 2 ^ 1 —x 1 cos -1 c— 2 x log (c—x) 

+2x log (1 —cx— V 1 — x 2 V 1 — c 2 ) 

=2(Mlog(t^g^) 

— 2 ^ 1 —x 1 cos _I c 

EVALUATION OF 
v>i= J {log[(x— x0 2 +(y— y,) 2 ] 

— log[(x— x^H (y+yi) 2 ]} (xi— c)dx. 


(fr-cf 


{log[(x-x,) 2 +(y-y,) 2 ] 


-log[(x-x,) 2 +(y+y.) 2 )}]; 

p fa — c) 2 dxi _ p (xi—c) 2 dxi _ _ f (cos 0— c) 2 d0 
Jc i/i(x — *i) ~ Jc l—x i l (x—x 1 )~ J X-COS0 

Xi = cos 6, y x = sin 9, dxj — — sin 0d0 

rvss? =sin9+(i " 2c)# " (a5 “ c)s i? de 

p d0 p d(7r+g) 

J c x— COS d Jc 37 -j - COS (tt + 0) 


COS 0 


1 . 1 —a: cos 0— v 1 — x 2 sin 

Vl“ log 


XT " X COS 0 j cos 0=c 

l— cx— y 1 — x 2 V 1 — c 2 


L_n i — x i i 

= Vi^L log ^T _log " 


Tj^logCl-cx-Vl-^Vl-c 2 ) 


=5 log (x-c) 


V 1 — X 2 

~-J ( p x = Vl— X 2 ^ — Vl — C 2 — (j? — 2 c) COS" 1 ! 
(x-cV 




(x— c)\ , vi 

-VT^ l08(x_<:) J 


— <p x = — Vl — C 2 ^l— X 1 — COS 'c (3; — 2c) V 1 — X 2 


X 


+ ( 37 — c) 2 log (1— C37— Vl — ^Vl— C 2 ) 

— (x— c) 2 log (x — c) 

EVALUATION OF T 3 
1 2tt <Px(% c)dx— -y/l c 2 f( a ._ c) ^f3pdx 

C € J 

cos _1 cJ* (x— c) (x— 2c) yl— x 2 dx 


+ ^ C -- log (1 — cx— V 1 —x 2 V 1 ~ c 2 ) 


— ^ J* (X— c) 8 dx — V 


r=?(x-c) 3 


dx 


4 V 1 — x 2 

— (x—c) 3 log (x— c)dx; x=cos 9, dx=— sin 

Y j[Vz(x— c)dx= Vl— c 2 J* ( cos sin ' 2 

4-cos -1 c J (cos 9— c) (cos 6— 2c) sin 2 0d0 
+ ■ aV ^ C L log (1 —cx— V 1 — x 2 Vl — c 2 ) 

-j J (x-c) 3 dx+ ^ ll ~ J (cos 0-c) 3 d0 

_ (x_c)_ log ( x — c ) J (aj— c ) 3 dx 
2 f J <Px(x — c)dx= — cos-'cj'cos 4 


Odd 


1 0d0 


-f^3c cos^c — Vl — c 2 + ^ ^ J* < 


cos 3 ede 


-|- ^cos _1 c — 2c 2 cos _1 c + c V 1 — c 2 — ^cVl — c 2 ^ J" cos 2 0d0 

+ (-3eco.^+Vi^+5^; 

-b^2c 2 cos _1 c — c V 1 — c 2 — — ^ |*d0 


cos 0d0 
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EVALUATION OF T 5 


, f <*>s 3 0 sin 0 , 3 ( 0 , sin 0 cos 6\] 
--cos c^ 4 '4 y 2 "‘ 2 )\ 

cos _1 c— | Vl-^^sin 0(cos 2 04-2) 


/ v 
4( cos -1 c — 2c 2 cos -1 c + - 


c V 1 — c 2 ^ ^ 0 + sin 0 cos 0 ^ 


2 * 


4^— 3c cos _1 c+ Vl — c 2 +^~^ — ~ 
-)-^2c 2 cos~*c — cVl — c 2 — ^ C ^0 

— COS c( g 7T+ 2~ 3?r j= — g 7T cos 'c 

J i 

^>r(x— c)dx 


2 ' 2 
sin I 


)„r 




rsjr-c’ . 3 cos _1 c . 3cv> — c ! ”l 

"L 4 ! 8 ' 8 J 

, ](^vT^?+ 2VT^) 


Vl — c‘ 


cos ] c— 2c 2 cos J cH — 


cyi — c 2 ^ ^COS~ 1 C-f c VI — C 2 ' 


XI 


1— cx— Vl— x 2 Vl — c 2 
x— c 


— 3ccos~ 1 c+ Vi ~~ c2 4- 3C -^Vl— c 2 

—(2 c 2 cos” *c — c V 1 — c 2 — ~~^4 — — ^ cos _1 c 

^os-'c^l— |+c 2 — 2c 2 J 

+ Vl— x 2 cos _, c^+^ C — c 3 — 2c—^-hc 3 4-3cH-c 

, c 3 c"l , c 2 (l — c 2 ) , (1— c 2 ) c 2 (l c 2 ) 

+ 4 8 j 4 + 2 8 

- (1 -c 2 ) — 3c2(1 4 ~ c2) = -f~ c 2 ) (cos” 1 c) 2 

_l_c Vl c 2 cos”*c (7 + 2 c 2 ) _( 1 ^!) ( 5 c 2 + 4 ) (=^ 3 ) 
4 o 


2(x— c) log : 

— 2Vl— x 2 cos” 1 c|dx=r 5 =— 2 J (x—c) log (x— c)dx 
+ 2 J (x—c) log (1 —cx— Vl — x 2 Vl —c 2 )dx 


—2 cos”‘c Vl— x 2 dx= — 2 ^ 0 ^ log — c) 

4- j* ( 3 — c)dx4-2 cos -I cJ*sir 
-f (x—c) 2 log (1— cx— Vl— x^l— c 2 ) 

^dx 


2 

sin 2 0d0 


-J(x-c)^- 


— c+x- 


V 1 — X 2 


Now 
(x-c ) 5 


- yi— x 2 Vi— c 2 


— c4-x 


Vl— c 2 

Vl —x 2 


1— cx— v 1— x 2 Vl— c 2 


J 

-/I- _ _ 

+ (i -^ 1 * - dI 

T h — — (x—c) 2 log (x— c)4-2 cos -1 cJ*sin 2 0d0 
4- (x—c) 2 log (1 — cx— Vl — X 2 Vl —C 2 ) 

+ Vl— ( cos c)d0 


2 cos -1 c 


(0— sin 0 cos 0) 4" Vl — £ 2 sin 0 


— cVl— c 2 0 


COS 0 = 1 
COS0 = C 

= — (1— c 2 ) — (cos -1 c) 2 +2cVl — c 2 cos~‘c 
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Positive directions of axes and angles (forces and moments) are shown by arrows 


Axis 

Designation 

Sym- 

bol 

Longitudinal 

Lateral _ 

X 

Y 

Z 

Normal. _ 



Force 
(parallel 
to axis) 
symbol 


X 

Y 

Z 


Moment about axis 


Designation 


Rolling 

Pitching. T __ 
Yawing 


Sym- 

bol 


L 

M 

N 


Positive 

direction 


Y- 

Z- 

X— 


-*z 

->x 

+Y 


Angle 


Designa- 

tion 


Roll 

Pitch 

Yaw 


Sym- 

bol 


Velocities 

Linear 
(compo- 
nent along 
axis) 

Angular 

u 

V 

V 

q 

w 

r 


Absolute coefficients of moment 

c-± n —— o-^- 

L ‘~qbS m ~qcS qbS 

(rolling) (pitching) (yawing) 


Angle of set of control surface (relative to neutral 
position), 5. (Indicate surface by proper subscript.) 


4. PROPELLER SYMBOLS 


D 

V 

p/D 

V' 

V, 

T 

Q 


Diameter 
Geometric pitch 
Pitch ratio 
Inflow velocity 
Slipstream velocity 


Thrust, absolute coefficient C T = — 

’ pn 2 D 4 

Torque, absolute coefficient 


P 

C. 

v 

n 


Power, absolute coefficient C P = 


pri^D 5 


pV 6 

Pn 2 


Speed-power coefficient=^- 
Efficiency 

Revolutions per second, rps 
Effective helix angle =tan -1 ^^^ 


5. NUMERICAL RELATIONS 


1 hp=76.04 kg-m/s=550 ft-lb/sec 
1 metric horsepower =0.9863 hp 
1 mph= 0.4470 mps 
1 mps =2.2369 mph 


1 lb=0.4536 kg 
1 kg =2. 2046 lb 
1 mi= 1,609.35 m=5,280 ft 
1 m=3.2808 ft 


U. S. GOVERNMENT PRINTING OFFICE! 1941 


51 


APPENDIX B 

MATLAB® PROGRAM FOR RE-COMPUTING RESULTS IN NACA 496 

This appendix contains a Matlab® program that implements the two-degree-of-freedom flutter 
equations and the two-degree-of-freedom solution method found in Appendix I of NACA 496. By 
executing this program, results in figures 5 through 11 and 15 through 17, found in Appendix II of NACA 
496, may be re-computed directly. The re-computation of results in figures 12 through 14 is discussed 
below. 

The program is structured as depicted in the simplified flow chart, below, on the left. There are four 
"groups" of re-computed figures that may be chosen, according to the table, below, on the right. The 
user inputs an integer from 1 to 4, thereby specifying a group. Then, based on the integer, the program 
selects a set of quantities [k, Xa, xp, rp, etc.), computes the coefficients of the equations of motion, 
solves the equations of motion and plots the results. The "cases" in the table are consistent with 
subsection titles in section VI of the main body of the present paper. 



The re-computation of results in figures 12 through 14 requires multiple executions of the program, with 
the user selecting group number 2 and supplying appropriately modified sets of quantities (see figure 
legends of figs. 12 through 14) for this group. 

During the normal execution of the solution method, the computed quantity F and the computed /7 s 
may be, at times, real and, at other times, complex. Only real values of F and the /7s represent physical 
solutions; complex values of F and the /7s represent non-physical solutions. 

The Matlab® program in this appendix produces both physical and non-physical solutions. When 
plotting, Matlab® does not distinguish between real quantities and complex quantities and, if a quantity 
is complex, Matlab® will plot its real part and issue a warning message in the command window. The 
user is cautioned to check the command window for such messages, and if present, carefully examine 
the variable F and the /7s to discern which values in the solution are real and which are complex. By so 
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doing the user may determine if a particular locus of points in a figure represents a physical solution or a 
non-physical solution. 


%% This program reproduces figs 5 --> 11 and 15 --> 17 from NACA 496 
% Figs 12 --> 14 require multiple executions of group number 2 with 
% appropriately modified input quantities 

o, 

o 

% Implemented using NACA 496 equations and 2D0F solution method 
% Boyd Perry, III 

% NASA-Langley Research Center 

% May 2015 

clc 
clear 
close all 

%% Define k and kinv 

N = 10001; 
oneoverkmax = 100; 

deltaoneoverk = oneoverkmax/ (N-l) ; 
for i=l : N 

kinv(i) = ( i-1 ) ^deltaoneoverk; 
k (i) = 1/kinv (i) ; 

end 

%% Define Quantities 
a = [0 -0.2 -0.4] ; 
c = [0 0.5] ; 

% For all calculations a = -0.4 and c = 0.5 

group = input ('Enter group number (1, 2 , 3, or 4) ') 

% Group number may be 1 , 2 , 3, or 4 

% Group number 1 is for the Standard Case and reproduces figs 5 --> 10 

% Group number 2 is for the Standard Case and reproduces fig 11 

% Group number 3 is for Wing A and reproduces fig 15 

% Group number 4 is for Wing B and reproduces figs 16 & 17 

% Input quantities for group number 1 

if group == 1 

kappa = 0.1; 

xa = 0.2; 

ra2 = 0.25; 


xb = 

: 1/80; 

o 

o 

Corresponds 

to 

curve 

(a) 

in 

figs 

7 

& 

8 

%xb 

= 1/40; 

o 

o 

Corresponds 

to 

curve 

(b) 

in 

figs 

7 

& 

8 

%xb 

= 1/160; 

o 

o 

Corresponds 

to 

curve 

(c) 

in 

figs 

7 

& 

8 


rb2 = 1/160; 
end 
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% Input quantitie for group number 2 

if group == 2 

kappa = 1/400; 

xa = 0.2; 

ra2 = 0.25; 

xb = 1/80; 

rb2 = 1/160; 

end 

% Input quantitie for group number 3 
if group == 3 
kappa = 1/416; 

%xa = 0.31; 
xa = 0.173; 

%xa = 0.038; 
ra2 = 0.33; 
xb = 1/80; 
rb2 = 1/160; 
omegaa = 7*2*pi; 
b = 2.5/12; 
end 

% Input quantitie for group number 4 

if group == 4 

kappa = 1/100; 

xa = 0.192; 

ra2 = 0.178; 

xb = 0.019; 

%xb = 0.01; 
rb2 = 0.0079; 
omegaa = 17.6*2*pi; 
omegah = 5.8*2*pi; 
omegab = 4.4^2^pi; 
b = 2.5/12; 
end 

% Other relationships 

ra = sqrt ( ra2 ) ; 
rb = sqrt ( rb2 ) ; 

%% Define array indices 

o 

o 

% For Ti constants : 

% Ti(kk) 

% kk ~ variation of quantity c 

o 

o 

% For R' arrays: 

% R’ ( j , kk) 

% j ~ variation of quantity a 
% kk ~ variation of quantity c 

o 

o 

% For F and G, real and imag parts of Theodorsen Function 
% F (i) and G (i) 
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i ~ variation of reduced frequency 


o 

o 

o 

o 

o 

o 

o, 

o 

o 

o 

o 

o 

o 

o 

o 

o 


For R" and I arrays : 

R M (i,j,kk) and I(i,j,kk) 
i ~ variation of reduced frequency 
j ~ variation of quantity a 
kk ~ variation of quantity c 


%% Define Ti 

for kk=l:2 % loop on quantity c [0, 0.5] 

c2 (kk) = c (kk) *c (kk) ; 

TiTerml = 1 + c2 (kk) ; 

TiTerm2 = 1 - c2 (kk) ; 

TiTerm3 = sqrt (TiTerm2 ) ; 

TiTerm4 = acos(c(kk)); 

TiTerm5 = 1/8 + c2 (kk) ; 

TiTerm6 = 7 + 2*c2 (kk) ; 

Tl (kk) = - ( 1/3 ) *TiTerm3* (2+c2 ( kk) ) + c ( kk) *TiTerm4 ; 

T2 (kk) = c ( kk) *TiTerm2 - TiTerm3*TiTerml*TiTerm4 + c ( kk) *TiTerm4 *TiTerm4 ; 
T3(kk) = -TiTerm5*TiTerm4*TiTerm4 + 

(1/4) * c (kk) *TiTerm3*acos (c (kk) ) *TiTerm6 - (1/8) *TiTerm2* (5*c2 (kk) +4) ; 

T4 (kk) = -TiTerm4 + c ( kk) *TiTerm3 ; 

T5 (kk) = -TiTerm2-TiTerm4 ^TiTerm4 + 2 ( kk) *TiTerm3*TiTerm4 ; 

T6 (kk) = T2 (kk) ; 

T7 (kk) = -TiTerm5*TiTerm4 + ( 1/8 ) ( kk) *TiTerm3*TiTerm6 ; 

T8 (kk) = - (1/3) *TiTerm3* (2^c2 (kk) +1) + c ( kk) *TiTerm4 ; 

T10 (kk) = TiTerm3 + TiTerm4; 

Til (kk) = TiTerm4 * ( 1-2 ( kk) ) + TiTerm3* (2-c (kk) ) ; 

T12 (kk) = TiTerm3* (2+c (kk) ) - acos (c (kk) ) * (2*c (kk) +1) ; 

end 

%% Define F and G 
for i=l : N 

JO (i) = besselj (0 , k (i) ) ; 

Jl(i) = besselj (1, k (i) ) ; 

Y0 (i) = bessely (0 , k (i) ) ; 

Yl (i) = bessely ( 1 , k (i ) ) ; 

TheoTerml (i) = J1 (i) + Y0 (i) ; 

TheoTerm2 (i) = Yl (i) - JO (i) ; 

Fnumer(i) = J1 (i) *TheoTerml (i) + Yl (i) *TheoTerm2 (i) ; 

Gnumer(i) = Yl (i) *Y0 (i) + Jl (i) *J0 (i) ; 

denom(i) = TheoTerml (i) A 2 + TheoTerm2 ( i ) A 2 ; 

F(i) = Fnumer ( i ) /denom ( i ) ; 

G(i) = -Gnumer ( i ) /denom (i); 


end 

%% Define R' Terms 

for kk=l : 2 % loop on quantity c [0, 0.5] 

for j =1 : 3 % loop on quantity a [0 -0.2 -0.4] 
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Raap(j,kk) = -ra2/kappa - (1/8 + a(j)*a(j)); 

Rabp(j,kk) = -rb2/kappa - ( c ( kk) -a ( j ) ) *xb/kappa + T7 (kk) /pi + (c(kk)- 

a ( j ) ) *T1 (kk) /pi; 

Rahp(j,kk) = -xa/kappa + a(j); 


Rbap ( j , kk) 
Rbbp ( j , kk) 
Rbhp ( j , kk) 


Rabp ( j , kk) ; 

-rb2/kappa + T3 (kk) / (pi*pi) ; 
-xb/kappa + Tl (kk) /pi; 


Reap ( j f kk) 
Rebp ( j , kk) 
Rchp ( j , kk) 


Rahp ( j f kk) ; 
Rbhp ( j , kk) ; 
-1/ kappa - 1; 


end 

end 


%% Define R" and I Terms 

% Equations (1) — > (9) and (11) — > (19) on p. 14 


for kk=l:2 % loop on quantity c [0, 0.5] 

p = - (1/3) * (1-c (kk) *c (kk) ) A (3/2) ; 

for j =1 : 3 % loop on quantity a [0 -0.2 -0.4] 

for i=l:N % loop on reduced frequency 


RITerml = (a(j)+0.5) ; 

RITerm2 = (0 . 5-a ( j ) ) *G (i) - kinv ( i ) *F ( i ) ; 

RITerm3 = Til (kk) *G (i) -2* kinv (i) *T10 (kk) *F (i) ; 


Raapp ( i , j , kk) = kinv (i) *2*RITerml*RITerm2 ; 

Rabpp (£, j , kk) = 

kinv (i) * (1/pi) * ( (T4 (kk) +T10 (kk) ) *kinv (i) +RITerml*RITerm3 ) ; 

Rahpp ( i , j , kk) = kinv (i) *2*RITerml*G (i) ; 

Rbapp ( i , j , kk) = -kinv(i)*(T12(kk)/pi)*RITerm2; 

Rbbpp ( i , j , kk) = -kinv (i) * (1/pi) * (1/pi) * ( (T12 (kk) / 2) *RITerm3- 
kinv(i) * ( T5 (kk) -T4 (kk) *T10 (kk) ) ) ; 

Rbhpp (i, j , kk) = -kinv(i)*(T12(kk)/pi)*G(i); 

Rcapp ( i r j r kk) = -kinv(i)^2^RITerm2; 

Rcbpp ( i r j r kk) = -kinv(i)*(l/pi)*RITerm3; 

Rchpp ( i , j , kk) = -kinv ( i ) *2 *G ( i ) ; 

RITerm4 = ( 0 . 5-a ( j ) ) *F ( i ) + kinv (i) *G (i) ; 

RITerm5 = Til (kk) *F (i) +2^ kinv (i) ^T10 (kk) *G (i) ; 


Iaa (i f j , kk) 
lab (i f j , kk) 
a ( j ) ) *T4 (kk) ) ; 

Iah (i, j , kk) 


-2 ^RITerml ^RITerm4 + ( 0 . 5— a ( j ) ) ; 

- ( 1/pi) * (RITerml*RITerm5 + 2*p + (0.5— 

-2^RITerml^F (i) ; 


Iba (i f j , kk) 
Ibb (i , j , kk) 
Ibh (i , j , kk) 


(T12 (kk) /pi) *RITerm4 + (1/pi) * (p-Tl (kk) -0.5^T4 (kk) ) 
(0.5/ (pi*pi) ) * ( T12 (kk) *RITerm5-T4 (kk) *T11 (kk) ) ; 

( T12 (kk) /pi) *F(i) ; 
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Ica (i, j , kk) = 
Icb ( i , j , kk) = 
Ich ( i , j , kk) = 

end 

end 

end 

%% Create R's and I's for 
j = 3; 


kk = 2; 


for i=l : N 


R1 (i) 
R2 (i) 
R3 (i) 

= Raap(j,kk) + 
= Rabp ( j , kk) + 
= Rahp ( j , kk) + 

R4 (i) 
R5 (i) 
R6 (i) 

= Rbap ( j , kk) + 
= Rbbp ( j , kk) + 
= Rbhp ( j , kk) + 

R7 (i) 
R8 (i) 
R9 (i) 

= Reap ( j , kk) + 
= Rcbp(j,kk) + 
= Rchp(j,kk) + 

111 (i) 
112 (i) 
113 (i) 

= Iaa (i, j , kk) 
= lab (i, j , kk) 
= Iah ( i , j , kk) 

114 (i) 

115 (i) 

116 (i) 

= Iba (i, j , kk) 
= Ibb (i , j , kk) 
= Ibh (i , j , kk) 

117 (i) 
118 (i) 
119 (i) 

= Ica (i, j , kk) 
= Icb (i, j , kk) 
= Ich (i, j , kk) 


2*RITerm4 + 1; 

(1/pi) * (RITerm5-T4 (kk) ) ; 
2*F(i) ; 


j=3 (a=-0 . 4 ) and kk=2 (c=0.5) 


Raapp ( i , j , kk) ; 
Rabpp (i, j , kk) ; 
Rahpp (i, j , kk) ; 

Rbapp (i, j , kk) ; 
Rbbpp (i, j , kk) ; 
Rbhpp ( i , j , kk) ; 

Rcapp (i, j , kk) ; 
Rcbpp (i, j , kk) ; 
Rchpp (i, j , kk) ; 


end 


%% Solve Case 3 Example Using Equations from Appendix I (Figs 5 & 6) 
if group == 1 


% Compute quadratic coefficients 
for i=l : N 

MRch(i) = Rl(i)*R5(i) - R2(i)*R4(i) - kinv (i) *kinv (i) * (111 (i) *115 (i) - 
112 (i) *114 (i) ) ; 

Mlch(i) = Rl(i)*H5(i) - R2(i)*I14(i) + Ill(i)*R5(i) - 112 (i) *R4 (i) / 
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A ( i ) 
B ( i ) 
C ( i ) 


MRch ( i ) *115 (i) *115 (i) - 
-Mich (i) * (R2 (i) *114 (i) + 
MRch (i) *111 (i) *111 (i) - 


Mich (i) *R5 (i) *115 (i) ; 

112 (i) *R4 (i) ) + 2*MRch(i) *111 (i) *115 (i) 
Mich (i) *R1 (i) *111 (i) ; 


end 


% Solve quadratic equation 
for i=2 : N 

Quad = [A (i) B (i) C (i) ] ; 
r = roots (Quad) ; 

OmegaAlphal (i) = r(l); 

0megaAlpha2 (i) = r(2); 

XI (i) = -Mich (i) / (OmegaAlphal (i) *115 (i) +111 (i) ) ; 

X2(i) = -Mich (i) / (0megaAlpha2 (i) *115 (i) +111 (i) ) ; 

FI (i) = kinv (i) *kinv (i) /XI (i) ; % This is actually the square of F 
F2 (i) = kinv ( i ) *kinv ( i ) /X2 ( i ) ; % This is actually the square of F 

end 

figure ( 5 ) 

plot (kinv, OmegaAlphal, 'r.', kinv, 0megaAlpha2, 'b.'); 
grid 

axis ([010 160] ) 

xlabel ( ' 1/k ' , ' FontWeight ' , ' bold ' ) ; 
ylabel ( ' Omega-Alpha ' , ' FontWeight ' , 1 bold ' ) ; 

title ('Fig 5, Case 3, Standard Case ',' FontWeight ',' bold ') ; 


figure ( 6 ) 

plot (OmegaAlphal , FI, 'r.', OmegaAlpha2, F2, 'b.'); 

grid 

axis ( [0 180 0 22] ) 

xlabel ( ' Omega-Alpha ' , ' FontWeight ' , ' bold ' ) ; 
ylabel ( ' F 1 , 1 FontWeight ' , ' bold ' ) ; 

title ('Fig 6 Case 3, Standard Case ',' FontWeight ', 'bold'); 


end 


%% Solve Case 2 Example Using Equations from Appendix I (Figs 7 & 8) 


if group == 1 


% Compute quadratic coefficients 
for i=l : N 


MRaa ( i ) = R5 ( i ) *R9 ( i ) - R6 ( i ) *R8 ( i ) - 

116 (i) *118 (i) ) ; 

Mlaa(i) = R5 ( i ) * 1 1 9 ( i ) - R6 ( i ) * 1 1 8 ( i ) 


kinv (i) *kinv (i) * (115 (i) *119 (i) - 
+ 115 (i) *R9 (i) - 116 (i) *R8 (i) ; 


A ( i ) 
B ( i ) 
C(i) 


MRaa (i) *119 (i) *119 (i) 
-Mlaa (i) * (R6 (i) *118 (i 
MRaa (i) *115 (i) *115 (i) 


- Mlaa (i) *R9 (i) *119 (i) ; 

+ 116 (i) *R8 (i) ) + 2*MRaa (i) *115 (i) *119 (i) 

- Mlaa (i) *R5 (i) *115 (i) ; 
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end 


% Solve quadratic equation 
for i=2 : N 

Quad = [ A ( i ) B(i) C(i)]; 
r = roots (Quad) ; 

OmegaBetal (i) = r(l); 

0megaBeta2 (i) = r(2); 

XI (i) = -Mlaa (i) / (OmegaBetal (i) *119 (i) +115 (i) ) ; 

X2(i) = -Mlaa (i) / (0megaBeta2 (i) *119 (i) +115 (i) ) ; 

FI (i) = kinv (i) *kinv (i) /XI (i) ; % This is actually the square of F 
F2 (i) = kinv ( i ) *kinv ( i ) /X2 ( i ) ; % This is actually the square of F 

end 

figure ( 7 ) 

plot (kinv, OmegaBetal, f r.', kinv, 0megaBeta2, 'b.'); 
grid 

axis ( [0 4 -0.008 0.014] ) 

xlabel ( 1 1/k ' , 1 FontWeight ' , ' bold ' ) ; 

ylabel ( ' Omega-Beta 1 , ' FontWeight 1 , 1 bold ' ) ; 

title ('Fig 7, Case 2, Standard Case ',' FontWeight ',' bold ') ; 

annotation (' textbox String ',{[' xb = ', num2str (xb) ]},' FontWeight bold ' , 
' FontSize ' , 10 , . . . 

' BackgroundColor ' , [ 1 1 1]); 

figure ( 8 ) 

plot (OmegaBetal , FI, 'r.', OmegaBeta2, F2, 'b.'); 

grid 

axis ( [0 0.014 0 1.5] ) 

xlabel ( ' Omega-Beta ' , ' FontWeight ' , ' bold ' ) ; 
ylabel ( ' F ' , ' FontWeight ' , ' bold ' ) ; 

title ('Fig 8, Case 2, Standard Case ',' FontWeight ',' bold ') ; 

annotation (' textbox ',' String ',{[' xb = ', num2str (xb) ]},' FontWeight ',' bold ' , 
' FontSize ' , 10 , . . . 

' BackgroundColor ',[ 1 1 1]); 


end 

%% Solve Case 1 Example Using Equations from Appendix I (Figs 9 & 10) 
if group == 1 

% Compute quadratic coefficients 
for i=l : N 

MRbb(i) = R9(i)*Rl(i) - R7 (i) *R3 (i) - kinv (i) *kinv (i) * (119 (i) *111 (i) - 
117 (1) *113 (1) ) ; 

Mlbb(i) = R9 (i) *111 (i) - R7 (i) *113 (i) + I19(i)*Rl(i) - 117 (i) *R3 (i) ; 

A ( i ) = MRbb (i) *111 (i) *111 (i) - Mibb (i) *R1 (i) *111 (i) ; 
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B ( i ) = -Mibb (i) * (R7 (i) *113 (i) + 117 (i) *R3 (i) ) + 2*MRbb (i) *119 (i) *111 (i) 
C ( i ) = MRbb (i) *119 (i) *119 (i) - Mibb (i) *R9 (i) *119 (i) ; 


end 

% Solve quadratic equation 
for i=2 : N 


Quad = [A (i) B (i) C (i) ] ; 
r = roots (Quad) ; 

OmegaHl (i) = r(l); 

0megaH2 (i) = r(2); 

XI (i) = -Mibb (i) / (OmegaHl (i) *111 (i) +119 (i) ) ; 

X2(i) = -Mibb (i) / (OmegaH2 (i) *111 (i) +119 (i) ) ; 

FI (i) = kinv ( i ) *kinv ( i ) /XI ( i ) ; % This is actually the square of F 
F2 (i) = kinv ( i ) *kinv ( i ) /X2 ( i ) ; % This is actually the square of F 

end 

figure (9) 

plot (kinv, OmegaHl, f r.', kinv, 0megaH2, 'b.'); 
grid 

axis ( [0 10 -30 120] ) 

xlabel ( ' 1/k ' , ' FontWeight ' , ' bold ' ) ; 

ylabel ( ' Omega-H ' , ' FontWeight ' , ' bold ' ) ; 

title ('Fig 9, Case 1, Standard Case ',' FontWeight ',' bold ') ; 


figure ( 10 ) 

plot (OmegaHl , FI, 'r.', OmegaH2, F2, 'b.'); 

grid 

axis ( [0 22 0 1.5] ) 

xlabel ( ' Omega-H ' , ' FontWeight ' , ' bold ' ) ; 
ylabel ( 1 F 1 , ' FontWeight ' , ' bold ' ) ; 

title ('Fig 10, Case 1, Standard Case ',' FontWeight ',' bold ') ; 
end 

%% Solve Case 1 Example Using Equations from Appendix I (Fig 11) 
if group == 2 


% Compute quadratic coefficients 
for i=l : N 

MRbb(i) = R9(i)*Rl(i) - R7 (i) *R3 (i) - kinv (i) *kinv (i) * (119 (i) *111 (i) - 
117 (i) *113 (i) ) ; 

Mlbb(i) = R9 ( i ) * 1 1 1 ( i ) - R7 (i) *113 (i) + I19(i)*Rl(i) - 117 (i) *R3 (i) ; 

A ( i ) = MRbb (i) *111 (i) *111 (i) - Mibb (i) *R1 (i) *111 (i) ; 

B ( i ) = -Mibb (i) * (R7 (i) *113 (i) + 117 (i) *R3 (i) ) + 2*MRbb (i) *119 (i) *111 (i) 
C ( i ) = MRbb (i) *119 (i) *119 (i) - Mibb ( i ) *R9 ( i ) *11 9 ( i ) ; 
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end 


% Solve quadratic equation 
for i=2 : N 

Quad = [A (i) B (i) C (i) ] ; 
r = roots (Quad) ; 


OmegaHl 

(1) = r (1) / 





OmegaH2 

(i) = r (2) ; 





XI (i) = 

-Mibb (i) / (OmegaHl (i) 

*m < 

i) +H9 (i) ) ; 



X2(i) = 

-Mibb (i) / (OmegaH2 (i) 

*m < 

i) +119 (i) ) ; 



Fl(i) = 

kinv (i) /sqrt (XI (i) ) ; 

% F 

as defined on 

P P . 

14-15 

F2(i) = 

kinv ( i ) / sqrt (X2 ( i ) ) ; 

% F 

as defined on 

PP. 

14-15 

omlom21 

(i) = ra*sqrt (OmegaHl 

(i) ) ; 




omlom22 

(i) = ra*sqrt (OmegaH2 

(i) ) ; 





end 

figure (11) 

plot (omlom21, FI, ' r.', omlom22, F2 , ' b.'); 

grid 

axis ( [ 0 1.6666666667 0 1.75]) 

xlabel (' omega-1 / omega-2 ’ , ' FontWeight bold ' ) ; 
ylabel ( 1 F ' , ' FontWeight ' , ' bold ' ) ; 

title ('Fig 11, Case 1, Standard Case (but with kappa = 
1/4 00 ) ' , ' FontWeight ' , ' bold ' ) ; 


end 


%% Solve Case 1 Example Using Equations from Appendix I (Fig 15) 
if group == 3 


% Compute quadratic coefficients 
for i=l : N 


MRbb (i) = R9 (i) 
117 (i) *113 (i) ) ; 
Mibb (i) = R9 (i) 

A ( i ) = MRbb ( i ) * 
B ( i ) = -MIbb(i) 
C ( i ) = MRbb ( i ) * 


*R1 ( i ) - R7 i)*R3 i - kinv(i)* 


*111 (i) - R7 (i) 

111 (i) *111 (i) - 
* (R7 (i) *113 (i) + 
119 (i) *119 (i) - 


113 (i) + 119 (i) 

Mibb (i) *R1 (i) *1 
117 (i) *R3 (i) ) + 
Mibb (i) *R9 (i) *1 


kinv ( i ) * 

(119 

d) *m <i) - 

*R1 (i) - 

117 

(i) *R3 (i) ; 

11 (i) ; 




2 *MRbb ( i ) *119 (i) *111 (i) 
19(1) ; 


end 


% Solve quadratic equation 
for 1=2 : N 


Quad = [ A ( 1 ) B(i) C(i)]; 
r = roots (Quad) ; 


61 


OmegaHl (i) = r(l); 

OmegaH2 (i) = r(2); 

XI (i) = -Mibb (i) / (OmegaHl (i) *111 (i) +119 (i) ) ; 

X2(i) = -Mibb (i) / (OmegaH2 (i) *111 (i) +119 (i) ) ; 

FI (i) = kinv (i) /sqrt (XI (i) ) ; % F as defined on pp . 14-15 
F2 (i) = kinv (i) /sqrt (X2 (i) ) ; % F as defined on pp . 14-15 

facl = b*omegaa*ra/sqrt ( kappa) ; 

fac2 = 12/39.37; % convert from fps to mps 

veil (i) = FI (i) *facl*fac2; 
vel2 (i) = F2 (i) *facl*fac2; 

omlom21 (i) = ra*sqrt (OmegaHl ( i )) ; 
omlom22 (i) = ra*sqrt (OmegaH2 ( i ) ) ; 


end 

figure ( 15 ) 

plot (omlom21 , veil, f r.', omlom22, vel2, 'b.'); 
grid 

axis ( [0 1.5 0 50] ) 

xlabel ( ' omega-h / omega-alpha FontWeight bold ' ) ; 
ylabel (' Velocity, meters per sec FontWeight bold ' ) ; 
title (' Fig 15, Case 1, Wing A FontWeight bold ') ; 

annotation (' textbox String ',{[' xa = ', num2str (xa) ]},' FontWeight bold ' , 
’ FontSize 1 , 10 , . . . 

1 BackgroundColor ' , [ 1 1 1]); 


end 

%% Solve Case 2 Example Using Equations from Appendix I (Fig 16) 
if group == 4 

% Compute quadratic coefficients 
for i=l : N 


MRaa(i) = R5(i)*R9(i) - R6(i)*R8(i) - kinv (i) *kinv (i) * (115 (i) *119 (i) - 
116 (i) *118 (i) ) ; 

Mlaa(i) = R5(i)*I19(i) - R6(i)*I18(i) + I15(i)*R9(i) - 116 (i) *R8 (i) ; 

A ( i ) = MRaa (i) *119 (i) *119 (i) - Mlaa (i) *R9 (i) *119 (i) ; 

B ( i ) = -Mlaa (i) * (R6 (i) *118 (i) + I 1 6 ( i ) *R8 ( i ) ) + 2*MRaa (i) *115 (i) *119 (i) ; 
C ( i ) = MRaa (i) *115 (i) *115 (i) - Mlaa (i) *R5 (i) *115 (i); 


end 


% Solve quadratic equation 
for i=2 : N 


Quad = [A (i) B (i) C (i) ] ; 
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r = roots (Quad) ; 

OmegaBetal (i) = r(l); 

OmegaBeta2 (i) = r(2); 

XI (i) = -Mlaa (i) / (OmegaBetal (i) *119 (i) +115 (i) ) ; 

X2 (i) = -Mlaa (i) / (OmegaBeta2 (i) *119 (i) +115 (i) ) ; 

FI (i) = kinv (i) /sqrt (XI (i) ) ; % F as defined on pp . 14-15 
F2 (i) = kinv (i) /sqrt (X2 (i) ) ; % F as defined on pp . 14-15 

facl = b*omegah/sqrt ( kappa) ; 

fac2 = 12/39.37; % convert from fps to mps 

veil (i) = FI (i) *facl*fac2; 
vel2 (i) = F2 (i) *facl*fac2; 

omlom21 (i) = sqrt (OmegaBetal (i) ) /rb; 
omlom22 (i) = sqrt (OmegaBeta2 (i) ) /rb; 


end 

figure (16) 

plot (omlom21 , veil, 'r.', omlom22, vel2, 'b.'); 
grid 

axis ( [0 1.8 0 50] ) 

xlabel (' omega-beta / omega-h FontWeight bold ') ; 
ylabel (' Velocity, meters per sec FontWeight bold ' ) ; 
title (' Fig 16, Case 2, Wing B 1 , 1 FontWeight * , ! bold ! ) ; 

annotation (' textbox String ',{[' xb = ', num2str (xb) ]},' FontWeight bold ' , 
’ FontSize 1 , 10 , . . . 

' BackgroundColor ’ , [ 1 1 1]); 


end 


%% Solve Case 3 Example Using Equations from Appendix I (Fig 17) 
if group == 4 

% Compute quadratic coefficients 


for i=l : N 

MRch(i) = Rl(i)*R5(i) - R2(i)*R4(i) - kinv (i) *kinv (i) * (111 (i) *115 (i) - 
112 (i) *114 (i) ) ; 

Mlch(i) = R1 (i) *115 (i) - R2 (i) *114 (i) + 111 (i) *R5 (i) - 112 (i) *R4 (i) ; 

A ( i ) = MRch(i) *115 (i) *115 (i) - Mich (i) *R5 (i) *115 (i) ; 

B ( i ) = -MIch(i) * (R2 (i) *114 (i) + 112 (i) *R4 (i) ) + 2*MRch (i) *111 (i) *115 (i) ; 
C(i) = MRch(i) *111 (i) *111 (i) - Mich (i) *R1 (i) *111 (i); 


end 


% Solve quadratic equation 
for i=2 : N 
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Quad = [ A ( i ) B(i) C(i)]; 
r = roots (Quad) ; 

OmegaAlphal (i) = r(l); 

OmegaAlpha2 (i) = r(2); 

XI (i) = -Mich (i) / (OmegaAlphal (i) *115 (i) +111 (i) ) ; 

X2 (i) = -Mich (i) / (OmegaAlpha2 (i) *115 (i) +111 (i) ) ; 

FI (i) = kinv ( i ) /sqrt (XI ( i ) ) ; % F as defined on pp . 14-15 
F2 (i) = kinv (i) /sqrt (X2 (i) ) ; % F as defined on pp . 14-15 

facl = b*omegab*rb/sqrt ( kappa) ; 

fac2 = 12/39.37; % convert from fps to mps 

veil (i) = FI (i) *facl*fac2; 
vel2 (i) = F2 (i) *facl*fac2; 

omlom21 (i) = sqrt (OmegaAlphal ( i ))*( rb/ra) ; 
omlom22 (i) = sqrt (OmegaAlpha2 ( i ))* (rb/ra); 


end 

figure ( 17 ) 

plot (omlom21 , veil, 'r.', omlom22, vel2, 'b.'); 
grid 

axis ( [0 2.8 0 40] ) 

xlabel (' omega-alpha / omega-beta FontWeight bold ') ; 
ylabel (' Velocity, meters per sec FontWeight bold ') ; 
title (' Fig 17, Case 3, Wing B FontWeight bold ') ; 

end 
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APPENDIX C 

ALTERNATE TWO-DEGREE-OF-FREEDOM SOLUTION METHODS 

Three alternate solution methods are offered, all of which are based on the two-degree-of-freedom 
flutter equations of NACA 496. Compared to the number of computations involved in the NACA 496 
solution method, these alternate solution methods are much more computationally intensive and would 
have been excluded from consideration in 1934 for that reason. But with today's computational 
resources, these alternate solution methods require only a few seconds of execution time on a desktop 
computer. The solution methods will be presented in order of increasing complexity. 

Although these alternate solution methods are appropriate for all three subcases mentioned in NACA 
496, the equations for Subcase 1, flexure - torsion, are chosen to illustrate the methods. The starting 
point for each alternate solution method is equation (XXIV) on page 12 of NACA 496, quadratic in X with 
complex coefficients, re-written here 

Cl h Cl a X 2 + (D. h A aa + D. a A ch )X + M b p = 0 (Cl) 


with companion equations 


X = 

£l h 

n a 


^r r o > r ^ 2 

[Cl) 

v vk ) 


("‘) 2 

\a) r r r y 

[CD 

a r a\ 2 

(C 4) 

x o) r r r ) 



where, for Subcase 1, the reference quantities, co r and ry, are chosen to be co a and r a , respectively. 

An illustrative example will be presented for each alternate solution method. The problem-specific 
quantities chosen for the example are those employed in NACA 496 for its "standard case": 


*•=0.1; c = 0.5; a = -0.4; x a = 0.2; r a 2 = 0.25; 


and with co a = 100, a>h - 50, and b- 1 

Using the appropriate values in equations (C3) and (C4), the values of Q, and f2 a are both 1. 

Because equation (Cl) is quadratic, zero, one, or two flutter solutions are possible. However, for the 
problem-specific quantities chosen, only one flutter solution is obtained. 

As expected and as will be seen, all three alternate solution methods will yield almost exactly the same 
value of flutter velocity, Vf, and almost exactly the same value of flutter reduced frequency, Ay, as those 
predicted by the two-degree-of-freedom solution method of NACA 496. The family of solutions for 
Subcase 1 for the standard case already exists in figures 9 and 10. When Q, = 1 is used in these figures 
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the problem-specific solution is obtained: i //= 173.26 feet per second (fps) and kf = 0.4355. The first 
row of Table Cl contains this result. 

On occasion throughout this appendix, the abbreviation ASM will be used to denote "alternate solution 
method." 


Alternate Solution Method No. 1 

This alternate solution method is directly analogous to the three-degree-of-freedom solution method 
discussed in NACA 496. It employs the computational shortcut of creating two equations, each with real 
coefficients, by separating the real and imaginary parts of equation (Cl), thereby eliminating complex 
arithmetic. The right-hand side of both new equations is zero. 

£l h £L a X 2 + (Q. h R aa + D. a R ch )X + M b p = 0 (C5) 

(Cl h I aa + Cl a I ch )X+ M' b p = 0 <C6) 

The first equation, (C5), is quadratic in X and is obtained from the real parts of A aa , A C h, and Mb# the 
second equation, (C6), is linear in X and is obtained from the imaginary parts. 

At this point in ASM No. 1, the artifice of treating X as a parameter, rather than as a known quantity, is 
employed. Equations (C5) and (C6) are each solved for their roots, identified herein asX Ri ,X R2 and X h 
for a large number of reduced frequencies. But, because from equation (C2), it is seen that X is 
proportional to the inverse square of velocity, only the real positive values of X Ri ,X Rz and Xj are 
retained. 


Next, the real positive roots are plotted as functions of the inverse of reduced frequency on the same 
set of axes. The intersection of either of the X R 's with the Xi identifies the value of X and the value of 
reduced frequency that simultaneously solve equations (C5) and (C6), and therefore also solves the 
original quadratic equation with complex coefficients, equation (Cl). These values are the flutter values 
-Xf and kf. 


With Xf and kf known, the artifice of treatingXas a parameter is abandoned and equation (C2) is 
employed, which when re-written to solve for flutter velocity, becomes 

_ _J_ 1 bcOgTa 

f JYf yfk kf 

Plugging Xf, kf, and the other quantities into equation (C7) yields the flutter velocity. 


(C7) 


Illustrative Example for Alternate Solution Method No. 1 . - ASM No. 1 was implemented in Matlab® 
where equations (C5) and (C6) were solved 10,000 times for values of the inverse of reduced frequency 
ranging from 0.001 to 10 in increments of 0.001. Figure Cl presents the results from this 
implementation: the ordinate represents the real positive roots of equations (C5) and (C6); the abscissa 
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is the inverse of reduced frequency. The solid black line represents the locus of the X Ri , the first real 
positive root of equation (C5); the dashed black line represents the locus of the X Rz , the second real 
positive root of equation (C5); and the solid red line represents the locus of the X If the real positive root 
of equation (C6). The circle represents the intersection of the X Ri locus with the Xj locus and was 
obtained using an interpolation scheme. 

By inverting the horizontal coordinate of the intersection, the flutter reduced frequency is obtained. By 
using it and the vertical coordinate in equation (C7), the flutter velocity is obtained. Both results are 
contained in the second row of Table Cl and are seen to agree extremely well with the NACA 496 results. 


Alternate Solution Method No. 2 

This alternate solution method is closely related to ASM No. 1 and also eliminates complex arithmetic 
via the same computational shortcut. But, rather than initially treating X as a parameter, as was done in 
ASM No. 1, in this solution method, X is assigned values and the flutter solution is obtained via 
systematic trial and error. 

ASM No. 2 employes an outer loop on velocity and an inner loop on reduced frequency and, using 
equation (C2), it computes a value of X at each passing. These X's are then substituted into the 
following equations 


+ ^a^ch) Z + ^b(S ~ Z 1 


(C8) 


(ftftJaa + ft a I ch) X + M bp ~ Z 2 


(C9) 


producing values of Zi and Z 2 at each passing. 

Recall that the flutter condition is obtained when Zi and Z 2 are both zero for the same values of v and k. 
However, unless luck intervenes or infinitesimally small increments in v and k are chosen for the outer 
and inner loops, it is almost assured that no combination of v and k will be found, which results in Zi and 
Z 2 being exactly zero simultaneously. This being the case, two different interpolation schemes are 
employed to find the combination of v and k that does result in Zi and Z 2 being exactly zero 
simultaneously. 

When each inner loop is complete (that is, for each velocity), curves of Zu and Z 2 as functions of reduced 
frequency are created. Because equation (C8) is quadratic, there can be zero, one, or two values of 
reduced frequency that cause Zi to be zero. Because equation (C9) is linear, there can be zero or one 
value of reduced frequency that causes Z 2 to be zero. 

The Zi and Z 2 curves may or may not pass through zero. For curves that do pass through zero, a first 
interpolation scheme is employed to identify the values of reduced frequency at the zero crossings and 
these values and their corresponding velocity are potential solutions and are stored in arrays. There can 
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be up to three arrays: two for the Z/s and one for the Z 2 s. These identified reduced frequencies are 
defined as k Zl ^,k Zl ^ and k Z2 , respectively. 


When the outer loop is complete, the reduced frequencies at the zero crossings, k z ,k z and k z , are 
plotted as functions of velocity. Each point on the k Zi ^snd k Zl ^ curves (that is, each reduced frequency 
and its companion velocity) corresponds to Zi being exactly equal to zero. Likewise, each point on the 
k Z2 curve corresponds to Z 2 being exactly equal to zero. The flutter condition is obtained when either 
the /c Zli curve or the k Zi ^ curve intersects the k Z2 curve. Via a second interpolation scheme, the values of 
v and k at the intersection are determined. These values are v/ and kf because they cause Zi and Z 2 to be 
zero simultaneously. 

Illustrative Example for Alternate Solution Method No. 2 . - This alternate solution method was also 
implemented in Matlab® where equations (C2), (C8), and (C9) were solved within nested loops. The 
inner loop varied reduced frequency from 0.01 to 2 in increments of 0.01; the outer loop varied velocity 
from 160 fps to 185 fps in increments of 0.5 fps. (The limits of the outer loop were chosen based on the 
knowledge of flutter velocity gained from the previous illustrative example.) 

Figure C2 presents a plot of Zu and Z 2 as functions of reduced frequency for a velocity of 160 fps, well 
below the flutter velocity. As can be seen Z lt the solid curve, intersects zero at two points, k Zl ^= 0.3706 
and k Zi ^= 0.5317, and Z 2 , the dashed curve, intersects zero at one point, k Z2 - 0.4740. These values of 
reduced frequency are combined with similar reduced frequencies for the other velocities to form the 
arrays that are plotted in the next figure. 

Figure C3 presents a plot of k z , k z , and , k z? as functions of velocity. As can be seen k z , the dashed 
curve, intersects the k z , the dashed-dotted curve, at the values of velocity and reduced frequency 
indicated in the third row of Table Cl. As with the results for ASM No. 1, those for ASM No. 2 are in 
excellent agreement with the results from NACA 496. 


Alternate Solution Method No. 3 

This alternate solution method is very straightforward, but uses complex arithmetic. It also initially 
employs the artifice of treating X as a parameter and later abandons it. 

Treating X as a parameter, equation (Cl), quadratic in X with complex coefficients, is solved for X directly 
for a large number of reduced frequencies. Because the coefficients are complex the resulting X's are 
also complex, but the X's are not complex conjugates. These complex X's are designated X and X 2 . 

With Xi and X 2 and their corresponding reduced frequencies known, the artifice is abandoned and 
equation (CIO) 


v = 


1 1 b(j 0 a r a 

y/xVx k 


(CIO) 
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is employed, yielding velocities, Vi correspondng to Xi, and V 2 , corresponding to X 2 , that are also 
complex and also not complex conjugates. 

All the complex X's, and therefore all the complex v's and their corresponding A's, satisfy equation (Cl) 
and are therfore all flutter solutions. However, because in actual practice the flutter velocity is real, not 
complex, the task now is to examine all of the complex v's to determine if any have a zero imaginary 
part, that is to determine if any are real. 

A convenient way to perform this examination is to plot the imaginary part of vj against the real part of 
Mi and the imaginary part of v 2 against the real part of m 2 , yielding two curves in the complex plain. Each 
point on each curve corresponds to a different value of reduced frequency. If one of the curves crosses 
the real axis, at the crossing, the imaginary part is zero, the real part is the flutter velocity, and the 
examination is complete. Interpolation is then employed, taking the complex velocities immediately 
before and immediately after the crossing to determine the value of the corresponding real part at the 
crossing, which is v/, and the corresponding value of reduced frequency at the crossing, which is Ay. 

Illustrative Example for Alternate Solution Method No. 3 . - This alternate solution method was also 
implemented in Matlab® where equations (Cl) and (CIO) were solved 10,000 times for reduced 
frequencies ranging from 0.005 to 50 in increments of 0.005. 

Figure C4 presents plots of the imaginary parts of vj and v 2 as functions of their respective real parts. 
The solid line is v 2 ; the dashed line, v 2 . By the nature of the curves, it is apparent that, as stated above, 
Mi and v 2 are not complex conjugates. As can be seen in the figure, the dashed curve crosses the real 
axis and the interpolated values of velocity and reduced frequency at the crossing are indicated in the 
fourth row of Table Cl. The results for ASM No. 3 are in perfect agreement with the results from NACA 
496. 


Table Cl 

Comparison of Flutter Predictions for 
Standard Case 


Solution Method 

Vf 

fps 

kf 

NACA 496 

173.26 

0.4355 

Alt. Soln. Meth. No. 1 

173.22 

0.4356 

Alt. Soln. Meth. No. 2 

173.21 

0.4358 

Alt. Soln. Meth. No. 3 

173.26 

0.4355 


69 


APPENDIX D 

THREE-DEGREE-OF-FREEDOM SOLUTION METHOD 


This appendix provides a solution to the three-degree-of-freedom flutter equation following the 
procedure recommended in NACA 496. Equation (XXI) on page 12 of NACA 496, re-written here as 
equation (Dl), is the three-degree-of-freedom flutter equation 

£l a £lpVL h X 3 + (Cl a Cl p A ch + (lp(l h A acc + £l h £l a A h p)X 2 

+(il a M acc + + £l h M ch )X + D = 0 


(Dl) 


with companion equations 


(br r a> r } 2 

(D2) 

\ vk ) 


(u a r a \ 2 
V co r r r ) 

(D3) 


(D4) 

\o) r r r J 


("*) 2 

\o) r r r y 

(D5) 


In equation (Dl) the A's, M's, and D are complex quantities, making equation (Dl) cubic in X with 
complex coefficients. 


The subsection entitled "Solution Methods" within section III of the main body of the present paper 
outlines the three-degree-of-freedom solution method. To recap, this solution method involves the 
artifice of treating X as a parameter and the computational shortcut of creating two equations, each 
with real coefficients, by separating the real and imaginary parts of equation (Dl), shown here as 
equations (D6) and (D7) 


tl a tlpil h X 3 + ( £l a ilpR ch + ilpil h R acc + a h £l a R h p)X 2 

+(n a ;iC + + n h M? h )x + d r = o 

{fl a £lpl ch + Uplift I aa + ^h^-ahp)X 2 

+ (A x M aa + ®.p M bp + ^h M ch) X + D l = 0 


(D6) 


(D7) 


Equations (D6) and (D7) are each solved for their roots, identified herein as X R ,X R ,X R ,X, and X, 2 , for 
a large number of reduced frequencies. But, because from equation (D2), it is seen that X is 
proportional to the inverse square of velocity, only the real positive values of X R ,X R ,X R ,X, and 
X/ 2 are of interest. 
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Next, all of the real positive roots are plotted as functions of the inverse of reduced frequency on the 
same set of axes. The intersection(s) of any of the X R 's with any of the X ! s identify(ies) the value(s) of X 
and the value(s) of reduced frequency that simultaneously solve equations (D6) and (D7), and therefore, 
also solves the original cubic equation with complex coefficients, equation (Dl). These values are the 
flutter values -Xf and kf. 

With Xf and the kf known, the artifice of treating X as a parameter is abandoned and equation (D2) is 
employed, which when re-written to solve for flutter velocity, becomes 

_ 1 1 bo) r v r (D8) 

Vf ~ 

Plugging Xf, kf, and the other quantities into equation (D8) yields the flutter velocity. 

The three-degree-of-freedom solution method is capable of finding zero, one, two, or three flutter 
modes, depending on the problem-specific quantities chosen. 


Illustrative Example for Three-Degree-of-Freedom Solution Method 

The problem-specific quantities chosen for this example are those employed in NACA 496 for its 
"standard case": 


*-=0.1; c = 0.5; a = -0.4; x a = 0.2; r a 2 = 0.25; xp=V S0 ; r/ = Yi 60 ; 

and with co a = 100; cop= 125; a>h = 50; b = 1; co r - 1; r r = 1 

This solution method was implemented in Matlab® where equations (D6) and (D7) were solved 10,000 
times for values of the inverse of reduced frequency ranging from 0.0005 to 5 in increments of 0.0005. 

Figure Dl presents the results from this implementation: the ordinate represents the real roots of 
equations (D6) and (D7); the abscissa is the inverse of reduced frequency. The three black lines (solid, 
dashed, dotted) represent the loci of the roots of equation (D6); the two red lines (solid, dashed) 
represent the loci of the roots of equation (D7); the open blue circle indicates the intersection of the 
locus of X R2 with the locus of X l± and is the flutter solution (a single flutter mode) for the quantities 
chosen. 

By inverting the horizontal coordinate of the intersection point, the flutter reduced frequency is 
obtained. By using it and the vertical coordinate of the intersection point in equation (D8), the flutter 
velocity is obtained. The results are contained in Table Dl. 
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Table Dl 

3D0F Flutter Prediction for 
Standard Case 


v f k f 

fps 


179.49 0.4476 


The problem-specific quantities for the illustrative examples in Appendices C and D of the present paper 
are the same, including the dimensional values of co a and cou- As stated in Appendix C of the present 
paper, alternate solution method no. 1 for the two-degree-of-freedom (2DOF) flutter equation is directly 
analogous to the three-degree-of-freedom (3DOF) solution method. Thus, in comparing these two 
illustrative examples, one might expect their results to be similar in some ways, which is, in fact, the 
case. 

Comparing first the values of flutter velocity and flutter reduced frequency from the two illustrative 
examples, from Tables Cl and Dl, it is seen that the respective values of Vf and Ay differ by only 3.5% and 
2.7%. These differences are attributed, obviously, to the presence of the third degree of freedom in the 
3DOF solution. 

Comparing next the loci of the roots in figures Cl and Dl, an obvious difference is the number of loci in 
each figure - three in the former and five in the latter. However, in figure Dl, if one ignores the loci of 
the third "real" root, X R3 , and the second "imaginary" root, X l2 , (those attributed to the presence of the 
third degree of freedom), the remaining loci agree very well in their general character with the loci in 
figure Cl. In each figure: (1) with increasing values of the inverse of reduced frequency, the loci of the 
first and second roots from the higher-order ("real") equation (shown in solid and dashed black) 
approach each other; (2) the locus of the root from the lower-order ("imaginary") equation (shown in 
solid red) lies between the other two loci; and (3) the locus of the "imaginary" root intersects the locus 
of the second "real" root. 


Results from 3DOF Solution Method Asymptotically Approach Results from 2DOF Solution Method 

Three examples are presented for the 3DOF solution method. In each example, two of the three modal 
frequencies (co a , cop, a>h) are fixed and the third is varied from its initial value to a very high value. In the 
limit, each example corresponds to one of the three subcases discussed in the main body of the present 
paper. For example, for Subcase 1, which has as degrees of freedom a and h, ry«and a>h are fixed while 
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cop is varied. With the value of the third modal frequency at least an order of magnitude higher than the 
highest value of the other two, the expectation is that the flutter velocities (v/) and flutter reduced 
frequencies (A/) obtained using the 3DOF solution method will asymptotically approach those quantities 
obtained using the appropriate 2DOF solution method. 

Results from 2DOF solution method . - The 2DOF solution methods are solved as described in the main 
body of this paper with quantities co r and r r defined according to which of the three 2DOF solution 
methods is employed. 

Each of these examples employs the problem-specific quantities from the NACA 496 "standard case": 

k=0.1; c = 0.5; a = -0.4; x a = 0.2; r a 2 = 0.25; xp=V 8 o; r/ = % iso 

The following table contains the natural frequencies and the corresponding values of frequency ratios 
for the three subcases. 


Table Dll 

Summary of Natural Frequencies and Frequency Ratios for Examples 



COa 

rps 

cop 

rps 

COh 

rps 

O a 

Op 

Qh 

Subcase 1 

100 

X 

50 

1 

X 

1 

Subcase 2 

X 

44.721 

50 

X 

0.005 

1 

Subcase 3 

100 

75 

X 

71.111 

1 

X 


To obtain the 2DOF flutter solutions for Subcase 1, figures 9 and 10 are used with Oh = 1; for Subcase 2, 
figures 7 and 8 are used with Op = 0.005; and for Subcase 3, figures 5 and 6 are used with O a = 71.111. 

Figures D2(a) and D2(b) contain close-up views of portions of figures 9 and 10; figures D3(a) and D3(b) 
contain close-up views of portions of figures 7 and 8; and figures D4(a) and D4(b) contain expanded 
views of figures 5 and 6. Following the convention established in the main body of the present paper, 
these curves are red, indicating they were created using present calculations. In each figure, the solid 
black line indicates the target value of its particular O. Intersections of the solid black line with the red 
curves produce the dashed black lines, which indicate the corresponding values of 1/A/ and F, from 
which Ay and v/, respectively, are obtained. (Recall that in figures 6, 8, and 10, although labelled "F" the 
quantity plotted is actually the square of the defined quantity, which has ramifications on the extraction 
of Vf from F.) 

In figures D2(a) and D2(b), there is only a single intersection of the black line with the red curves, 
indicating only one flutter mode is predicted for Subcase 1; but in figures D3(a) and D3(b) and in figures 
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D4(a) and D4(b), because there are two intersections of the black line with the red curves, two flutter 
modes each are predicted for Subcases 2 and 3. 

The following table presents a summary of the 2DOF flutter predictions: 


Table Dill 

2DOF Flutter Predictions for Standard Case 


Mode 1 Mode 2 



Vf 

fps 

kf 

Vf 

fps 

kf 

Subcase 1 

173.26 

0.4355 

X 

X 

Subcase 2 

19.521 

2.587 

120.65 

0.4727 

Subcase 3 

14.668 

8.045 

234.05 

0.4458 


Results from 3DOF solution method . - Figures D5, D6, and D7 contain the results of many 3DOF flutter 
solutions. These figures are comprised of semi-log plots of v/and kf as functions of the co appropriate to 
the subcase. The solid blue circles represent the values of v/ and kf predicted using the 3DOF solution 
method, where each corresponding pair came directly from an intersection such as the one indicated by 
the open blue circle in figure Dl. The dashed red lines in the figures represent the values of v/and kf 
predicted using the 2DOF solution method and contained in Table Dili. 

The nominal starting value for variable frequencies cop, co a , and coh in figures D5 through D7 was 100 rps. 
However, for Subcases 1 and 3, the fixed value of co a is also 100 rps. For these two subcases, the flutter 
modes that ultimately became the asymptotic solutions did not develop until their respective variable 
frequencies were above 100 rps: 125 rps for Subcase 1; 175 rps for Subcase 3. 

As can be seen in figure D5, the 3DOF solution method predicts a single flutter mode, as did the 2DOF 
solution method. Also, with increasing values of cop, the circles approach the dashed lines, indicating 
that, in the limit, the values of v/ and kf predicted by the 3DOF solution method approach those 
predicted by the 2DOF solution method, thereby confirming the initial expectation. 

As can be seen in figure D6, the 3DOF solution method predicts two flutter modes, as did the 2DOF 
solution method. In each plot, with increasing values of co a , the circles approach the dashed lines, again 
indicating that, in the limit, the values of Vf and /(/predicted by the 3DOF solution method approach 
those predicted by the 2DOF solution method, and again confirming the initial expectation. The 
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comparatively higher density of circles for this subcase between co a = 100 rps and co a = 250 rps was 
necessary in order to capture the region of high curvature and the inflection in figure D6(c). 

As can be seen in figure D7, the 3DOF solution method again predicts two flutter modes, as did the 
2DOF solution method. In each plot, with increasing values of a>h, the circles approach the dashed lines, 
once again indicating that, in the limit, the values of v/ and /(/predicted by the 3DOF solution method 
approach those predicted by the 2DOF solution method, and once again confirming the initial 
expectation. 

For each subcase, in the limit, the 3DOF solution method predicts the same number of flutter modes 
and the same values of flutter velocity and flutter reduced frequency as did the 2DOF solution method. 
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TABLE I. - TRIPPING POINTS IN NACA 496 


No. 

Tripping Point 

Identified herein 
on page - 

Occurs in 
NACA 496 - 

1 

Three versions of NACA 
496 dated 1934, 1940, 
and 1949. 

4 


2 

No errata for some of 
the errors in 1934 
version that were 
corrected in 1940 and 
1949 versions. 

4 


3 

Assignment of the 
same equation 
numbers, (1) - (9) and 
(11) - (19), to different 
expressions. 

12 

pp. 12, 14, 15 

4 

The use of the same 
symbol, /, y , to represent 
a quantity and that 
quantity divided by 
reduced frequency. 

13 

pp. 12 & 14 

5 

The use of the same 
symbol, F, to represent 
a quantity and the 
square of that quantity. 

16 

figs 6, 8, 10 - 
the square of 
the quantity 

figs 11-14- 
the quantity 

6 

Quantity D is never 
defined. 

18 

D is plotted 
in fig 12 

7 

Omission of values for 

coh and cop . 

23 

these values 
required to 
produce figs 
16 and 17 
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TABLE II. - IDENTIFIED ERRORS IN NACA 496 


No. 

Identified Error 

Identified herein 
on page(s) - 

Occurs in 
NACA 496 
on page - 

1 

The /?" h row of Table III corresponding to 
a = -0.4 and c = 0.5 is in error. 

15 

17 

2 

Expression^ = v R r~ is incorrect. 

-+a 

Should be v D = v R . 

19 

16 

3 

As presented in figure 13, the values 
comprising curve B are the square roots 
of the proper values. 

20 

20 

4 

Transcription errors: 

As presented in figure 14, correct values 
of F were paired with non-corresponding 
values of x a - that led to incorrect curves 
A, B, and C. 

22 

20 
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TABLE III. - COMPUTATIONAL SHORTCUTS IN NACA 496 


No. Computational Shortcut 

Benefit Identified herein Occurs in 

on page(s) - NACA 496 

on page(s) - 

1 Separating an equation with 
complex coefficients into its 
real and imaginary parts, 
yielding two equations with 
real coefficients. 

Avoids complex arithmetic. 9,10,66,67, 12 and 13 

and 70 

2 Eliminating the factor 1/b 
from equation (10/C), 
resulting in equations (12/3), 
(12/6), (12/9), (12/13), 
(12/16), and (12/19). 

Relieves engineers of extra, 12 12 

but unnecessary, 

multiplications in the hand- 

calculation of these 

expressions 

3 Eliminating the factor 1/k 
from equations (12/11) 
through (12/19), resulting in 
equations (14/11) through 
(14/19). 

Relieves engineers of extra, 12 14 

but unnecessary, 

multiplications in the hand- 

calculation of these 

expressions 
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Figure 1. - Parameters of the airfoil-aileron combination 
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c, location of aileron hinge c, location of aileron hinge c, location of aileron hinge 





- 1.0 - 0.5 0 0.5 1.0 - 1.0 - 0.5 0 0.5 1.0 

c, location of aileron hinge c, location of aileron hinge 


Figure 2. - Constants resulting from the integration of velocity potentials. 
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(a) Elements from pitching moment equation. 

Figure 4. - Reduced-frequency-dependent portions of the elements of the equations of motion. 
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(b) Elements from hinge moment equation. 

Figure 4. - Reduced-frequency-dependent portions of the elements of the equations of motion. 
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(c) Elements from vertical force equation. 

Figure 4. - Reduced-frequency-dependent portions of the elements of the equations of motion. 
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(d) Element f?"h re-drawn. 

Figure 4. - Reduced-frequency-dependent portions of the elements of the equations of motion. 
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Figure 5. - Case 3, Torsion-aileron (a, fi)\ Standard case. Showing Q a against V k . 





Figure 6. - Case 3, Torsion-aileron (a, ft): Standard case. Showing flutter factor 

F against 
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Figure 7. - Case 2, Aileron-deflection (J3, h ): (a) Standard case, (b), (c), (d) indicate 
dependency on Xp. Case (d), Xp = -0.004, reduces to a point. 



Figure 8. - Case 2, Aileron-deflection (J3, h): Final curves giving flutter factor F 
against Op corresponding to cases shown in figure 7. 
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Figure 9. - Case 1, Flexure-torsion ( h , a): Standard case. Showing Q h against V k 



Figure 10. - Case 1, Flexure-torsion ( h , a): Standard case. Showing flutter factor 

F against 
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CJ 2 /OJ z 


Figure 11 . - Case 1 , Flexure-torsion ( h , a): Showing dependency of F on <o h /(o a . The 
upper curve is experimental. Airfoil with r = Vi; a = -0.4; x = 0.2; 4k = 0.01; © j/© 2 variable. 
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Figure 12. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on location 
of axis of rotation a. Airfoil with r = !4; x = 0.2; k = %; co ( /(o 2 = V 6 ; a variable. 
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Figure 13. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on the radius 

of gyration r a = r. 

A, airfoil with a = -0.4; k=%;x = 0.2; ooj/oe^ = 34; r variable. 

B, airfoil with a = -0.4; k— %; x = 0.2; ooj/a^ = 1.00; r variable. 

Figure 13A. - Curve A 
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Figure 13. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on the radius 

of gyration r a = r. 

A, airfoil with a = -0.4; k='A;x = 0.2; co,/^ = %; r variable. 

B, airfoil with a = -0.4; k = %; x = 0.2; coj/o^ = 1.00; r variable. 


(a) Comparison of present calculation and Curve B of NACA 496 
Figure 13B. - Curve B 
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Figure 13. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on the radius 

of gyration r a = r. 

A, airfoil with a = -0.4; k='A;x = 0.2; co,/^ = %; r variable. 

B, airfoil with a = -0.4; k='A\x = 0.2; coj/o^ = 1.00; r variable. 


(b) Comparison of present calculation and square of Curve B of NACA 496 

Figure 13B. - Curve B 
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Figure 13. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on the radius 

of gyration r a = r. 

A, airfoil with a = -0.4; k='A;x = 0.2; co,/^ = %; r variable. 

B, airfoil with a = -0.4; k='A\x = 0.2; coj/o^ = 1.00; r variable. 


(c) Comparison of present calculation, independent calculation, and square of Curve B of NACA 496 

Figure 13B. - Curve B 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k = %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k = V 100 ; ©j l(o 2 = 1; x variable. 

(a) Comparison of present calculation and Curve A of NACA 496 
Figure 14A. - Curve A 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k= %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k= V 100 ; ©j l(o 2 = 1; x variable. 

(b) Comparison of present calculation and rotated Curve A of NACA 496 

Figure 14A. - Curve A 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = Vi; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k= %; c 0 j/g )2 = 'A; x variable. 

C, airfoil with r = Vi ; a = -0.4; k= V 100 ; ©j l(o 2 = 1; x variable. 

(c) Comparison of present calculation, independent calculation, and rotated Curve A of NACA 496 

Figure 14A. - Curve A 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k= %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k= V 100 ; ©j l(o 2 = 1; x variable. 

(a) Comparison of present calculation and Curve B of NACA 496 
Figure 14B. - Curve B 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k = %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k = V 100 ; ©j l(o 2 = 1; x variable. 

(b) Comparison of present calculation and rotated Curve B of NACA 496 

Figure 14B. - Curve B 


99 



Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = Vi; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k = %; C0j/g)2 = 'A', x variable. 

C, airfoil with r = Vi ; a = -0.4; k = V 100 ; ©j l(o 2 = 1; x variable. 

(c) Comparison of present calculation, independent calculation, and rotated Curve B of NACA 496 

Figure 14B. - Curve B 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; 0 O|/a) 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k = %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k = V 100 ; ©j l(o 2 = 1; x variable. 


(a) Comparison of present calculation and Curve C of NACA 496 
Figure 14C. - Curve C 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k = %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k = V 100 ; ©j l(o 2 = 1; x variable. 


(b) Comparison of present calculation and rotated Curve C of NACA 496 

Figure 14C. - Curve C 
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Figure 14. - Case 1, Flexure-torsion ( h , a): Showing dependency of F on x a , the 

location of the center of gravity. 

A, airfoil with r = !4; a = -0.4; k= V 400 ; (o,/co 2 = V 6 ; x variable. 

B, airfoil with r = 'A; a = -0.4; k = %; C0j/g)2 = 'A', * variable. 

C, airfoil with r = Vi ; a = -0.4; k = V 100 ; ©j l(o 2 = 1; x variable. 


(c) Comparison of present calculation, independent calculation, and rotated Curve C of NACA 496 

Figure 14C. - Curve C 
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Figure 15. - Case 1. Wing A. Theoretical and experimental curves giving flutter 
velocity v against frequency ratio oo h /ct) a . Deflection-torsion. 
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Figure 16. - Case 2. Wing B. Theoretical and experimental curves giving flutter 
velocity v against frequency ratio G)p/G) h . Aileron-deflection (/?, h). 
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Figure 17. - Case 3. Theoretical curve giving flutter velocity against the fre- 
quency ratio co a /cOp. The experimental unstable are is indefinite due to the im- 
portance of internal friction at very small velocities. Torsion-aileron (a, P). 
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Figure Cl - Illustrative example for Alternate Solution Method No. 1. 
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Figure C2 - Illustrative example for Alternate Solution Method No. 2 
-Z I and Z 2 as functions of reduced frequency at a velocity of 160 fps. 



Figure C3 - Illustrative example for Alternate Solution Method No. 2 
- Reduced frequencies as a function of velocity. 
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Figure D1 - Illustrative example for Three-Degree-of-Freedom Solution 
Method; "standard case" quantities chosen, with co a = 100, a>p- 125, a> h = 50. 
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Figure D2 - Results from 2DOF Solution Method, Subcase 1; 
co a = 100, co h - 50, r a 2 = 0.25, £2 h = 1. 




(a) Recreation of Figure 7 (b) Recreation of Figure 8 

Figure D3 - Results from 2DOF Solution Method, Subcase 2; 
a>p = 44.721, co h = 50, r fi 2 = 1/160, x fj = 1/80, Q p = 0.005. 




(a) Recreation of Figure 5 (b) Recreation of Figure 6 

Figure D4 - Results from 2DOF Solution Method, Subcase 3; 
co a = 100, (o p = 75, r a 2 = 0.25, r p 2 = 1/160, x p = 1/80, Q a = 71.111. 
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Figure D5 - Results from 3DOF Solution Method, "Subcase 1"; 
co a - 100, co h - 50, ^variable, "standard case". 
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Figure D6 - Results from 3DOF Solution Method, "Subcase 2"; 
<^=44.721, co h - 50, = 1/160, Xp = 1/80, co a variable. 
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(b) Mode 1 flutter reduced frequency vs co h 
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Figure D7 - Results from 3DOF Solution Method, "Subcase 3"; 
co a = 100, o)p= 75, r 2 = 0.25, r/ = 1/160, x R = 1/80, co h variable. 
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